AFIT/DS/ENP/96-01 


CHARACTERISTIC  SPATIAL  QUADRATURES 
FOR  DISCRETE  ORDINATES 
NEUTRAL  PARTICLE  TRANSPORT 
ON  ARBITRARY  TETRAHEDRAL  MESHES 

DISSERTATION 

Charles  R.  Brennan,  Captain,  USAF 
AFIT/DS/ENP/96-01 


W60718 118 


Approved  for  public  release;  distribution  unlimited 


AFIT/DS/ENP/96-01 


CHARACTERISTIC  SPATIAL  QUADRATURES 
FOR  DISCRETE  ORDINATES 
NEUTRAL  PARTICLE  TRANSPORT 
ON  ARBITRARY  TETRAHEDRAL  MESHES 

DISSERTATION 


Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 
Air  University  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 
Doctor  of  Philosophy 


Charles  R.  Brennan,  B.S.,  M  S. 
Captain,  USAF 
June  1996 


Approved  for  Public  release;  distribution  unlimited 


AFIT/DS/ENP/96-oi 


CHARACTERISTIC  SPATIAL  QUADRATURES  FOR  DISCRETE  ORDINATES 
NEUTRAL  PARTICLE  TRANSPORT  ON  ARBITRARY  TETRAHEDRAL  MESHES 


Approved: 


Charles  R  Brennan,  B.S.,  M.S. 
Captain,  USAF 


Kirk  A.  Mathews 

Associate  Professor  of  Nuclear  Engineering 
Research  Advisor 


L  C wlbrfi 


David  L.  Coulliette 

Assistant  Professor  of  Applied  Mathematics 


7 5  /U  y  U 


Todd  lA.Wareing 
Radiation  Transport  T  earn 
Los  Alamos  National  Laboratory 


William  P.  Baker 

Associate  Professor  of  Applied  Mathematics 
Dean's  Representative 


2  5  A1 


Dean,  Graduate  School  of  Engineering 


Acknowledgments 


I  am  eternally  grateful  to  the  many  people  who  assisted  me  during  my  research.  I 
owe  a  large  debt  to  Captains  Glenn  Sjonden,  Dennis  Miller  and  Bryan  Minor  for  building 
the  foundation  upon  which  my  work  rests.  I  wish  to  thank  Dr.  David  Weeks  for  all  of  his 
patience  when  I  bombarded  him  with  questions  in  the  early  stages  of  the  project.  I'm  also 
extremely  grateful  to  Dr.  Wally  Walters,  an  original  committee  member  who  passed  away 
during  my  research,  for  his  encouragement  and  help  in  defining  a  research  topic.  I  am 
indebted  to  the  members  of  my  committee.  Dr.  William  Bailey,  Lt.  Col.  Dave  Coulliette, 
and  Dr.  Todd  Wareing,  for  their  advice,  support,  and  encouragement.  I'm  especially 
grateful  to  my  research  advisor.  Dr.  Kirk  Mathews,  for  bravely  agreeing  to  take  me  on  as 
his  student,  for  laying  much  of  the  groundwork  leading  up  to  this  project,  for  his  advice 
and  assistance  in  defining  the  research  topic,  for  displaying  an  infinite  amount  of  patience, 
and  for  guiding  me  (successfully)  through  to  the  end.  Most  importantly,  I  thank  my  wife 
Fran,  and  my  children.  Tiffany  and  Meaghan.  None  of  this  would  have  been  possible 
without  their  patient  understanding  and  unending  support. 


Table  of  Contents 


Acknowledgments .  iii 

List  of  Figures . vii 

List  of  Tables .  xi 

Abstract .  xi 

I.  Introduction  . 1 

Background . 2 

Diffusion  Approximation  . 3 

Monte  Carlo  Methods  . 4 

Deterministic  Methods . 4 

Purpose . 9 

Scope . 10 

Organization  . 12 

II.  Derivation  of  the  Characteristic  Spatial  Quadratures . 13 

The  Boltzmann  Transport  Equation . 13 

Discrete  Ordinates  Formulation  . 16 

Boundary  Conditions  . 20 

Discrete  Ordinates  Solution . 21 

Role  of  the  Spatial  Quadrature . 24 

Moments  of  the  Boltzmann  Transport  Equation . 28 

Characteristic  Solutions  to  the  Boltzmann  Transport  Equation . 29 

Step  Characteristic . 29 

Linear  Characteristic . 30 

Exponential  Characteristic  . 31 

LC  and  EC  Variants  -  Constant  Face  Methods . 32 

Discretization  Errors  . 33 

III.  Derivation  of  the  Characteristic  Methods 

for  Arbitrary  Tetrahedral  Meshes . 35 

Application  of  the  Characteristic  Methods  to  Arbitrary  Tetrahedral  Meshes  ...  35 

Splitting  of  Arbitrary  Tetrahedra  into  Sub-Tetrahedra  . 36 

Definition  of  “Arbitrary”  Tetrahedra . 40 

Coordinate  Systems . 40 

Mathematical  Preliminaries:  The  Exponential  Moments  Functions  . 44 

Definition  and  Relationships . 45 


IV 


Derivation  of  the  Characteristic  Methods  for  Meshes  of  Arbitrary  Tetrahedra  .  .  47 

Step  Characteristic  Method  (SC)  . 47 

Linear  Characteristic  Method  (LC) . 53 

Exponential  Characteristic  Method  (EC) . 63 

Constant  Face  Variations  . 68 

Constant  Linear  Characteristic  . 68 

Constant  Exponential  Characteristic  . 69 

Balance  Equations . 70 

IV.  Arbitrary  Tetrahedra  Code  Algorithm . 74 

Initialization . 74 

Iteration . 74 

Pseudocode  Algorithm . 78 

V.  Testing  . 81 

Test  Problem  1  -  Unit  Source  Cube  with  Vacuum  Boundaries  . 83 

Convergence  . 85 

Mesh  Sensitivity . 96 

Execution  Speed . 106 

Test  Problem  2  -  Cube  Problem  with  Isotropic  Incident  Current 

and  Slab  Symmetry . 109 

Test  Problem  3  -  Cylinder-in-a-Cube,  Isotropic  Cylinder  Source . 117 

Test  Problem  4  -  Cylinder-in-a-Cube,  Isotropic  Incident  Current  . 125 

VI.  Summary/Recommendations . 135 

Summary  of  Results . 136 

Recommendations  for  Future  Efforts . 139 

Appendix  A.  Exponential  Moment  Functions  . 143 

One-Dimensional  Exponential  Moment  Functions  . 143 

Two-Dimensional  Exponential  Moment  Functions . 143 

Three-Dimensional  Exponential  Moment  Functions . 146 

Four-Dimensional  Exponential  Moment  Functions . 150 

Appendix  B.  Root  Solving  . 158 

Newton's  Method . 158 

Calculation  of  the  Input  Flux  Distribution  Coefficients  -  2D  System . 159 

Derivation  of  2D  System  Using  Multi-Dimensional 

Moment  Functions .  159 

Note  on  Notation . 160 


v 


Derivation  of  2D  System  Using  One-Dimensional 

Moments  Functions  . 162 

Initial  Guess  Equations  for  the  2D  System  . 164 

Computation  of  Distribution  Coefficients  from  Root-Solved  Values  .168 

Calculation  of  the  Source  Distribution  Coefficients:  3D  System . 171 

Derivation  of  3D  System  Using  Multi-Dimensional 

Moment  Functions .  171 

Derivation  of  3D  System  Using  One-Dimensional 

Moment  Functions .  173 

Initial  Guess  Equations  for  the  3D  System  . 174 

Computation  of  Distribution  Coefficients  from  Root  Solved  Values  .  .  179 

Appendix  C.  Tetrahedral  Splitting .  182 

Construction  of  the  Local  Sub-Tetrahedron  Coordinate  System . 182 

Case  4  Tetrahedra  Splitting  . 185 

Sub-Tetrahedron  1  (Configuration  1) . 188 

Sub-Tetrahedron  2  (Configuration  1) . 193 

Sub-Tetrahedron  3  (Configuration  1) . 194 

Sub-Tetrahedron  4  (Configuration  1) . 194 

Configuration  2  .  195 

Case  3  (-3)  Tetrahedra  Splitting . 195 

Case  2  (-2)  Tetrahedra  Splitting . 201 

Case  1  Sub-tetrahedra  . 202 

Appendix  D.  Conservation  Equations . 203 

Mesh  Generation . 206 

NASTRAN  Output  File . 207 

Sample  TetSN  Input  File  . 208 

Bibliography  . 209 

Vita . 214 


vi 


List  of  Figures 


Figure  1.  Problem  1  Tetrahedron  in  the  Local  uvw  Coordinate  System,  . 38 

Figure  2.  Tetrahedra  Splitting  and  Resulting  Number  of  Problem  1  Sub-Tetrahedra.  ..  39 

Figure3.  Orientation  of  Global  and  Centroid  Coordinate  Systems . 41 

Figure  4.  Global  Face  Coordinate  System . 42 

Figure  5.  Coordinate  System  Used  to  Compute  the  Coefficients  of  the 

Source  Distribution . 43 

Figure  6.  First  Three  Mesh  Refinements  Test  Problem  1  . . 84 

Figure  7.  Problem  Average  Flux  for  Test  Problem  1 . 87 

Figure  8.  Problem  Average  Flux  for  Test  Problem  1 . 

(Enlargement  of  fine  mesh  results) . 87 

Figure  9.  Exiting  Face  Flux  for  Test  Problem  1 . 88 

Figure  10.  Exiting  Face  Current  for  Test  Problem  1 . 

(Enlargement  of  Fine  Mesh  Results) . 88 

Figure  1 1 .  Relative  Error  in  Test  Problem  1  Problem  Average  Flux . 90 

Figure  12.  Relative  Error  in  Test  Problem  1  Exiting  Face  Current . 91 

Figure  13.  L-Infinity  Norm  (Maximum  Relative  Error)  of  Test  Problem  1 

Exiting  Face  Flux  Distribution . 93 

Figure  14.  L-l  Norm  (Average  Relative  Error)  of  Test  Problem  1 

Exiting  Face  Flux  Distribution . 93 

Figure  15.  Relative  Error  Between  TetSN  EC  and  MCNP  Solutions  for 

Test  Problem  1  Problem  Average  Flux . 95 

Figure  16.  Relative  Error  Between  TetSN  LC  and  MCNP  Solutions 

for  test  Problem  1  Problem  Average  Flux . 95 

vii 


Figure  17.  Relative  Error  Between  TetSN  SC  and  MCNP 

Solutions  for  Test  Problem  1  Problem  Average  Flux . 96 

Figure  18.  Exterior  View  of  25%  Variation  Meshes . 98 

Figure  19.  Exterior  View  of  the  50%  Variation  Meshes . 98 

Figure  20.  Cumulative  Distribution  Function  of  Tetrahedron  Volume 

Ratios  for  the  25%  Variation  Meshes . 99 

Figure  21 .  Cumulative  Distribution  Function  of  Tetrahedron  Volume 

Ratios  for  the  50%  Variation  Meshes . 100 

Figure  22.  Exterior  Mesh  View  of  Structured  Mesh,  Mesh  #4  and  Mesh  #5 

with  Worst  and  Best  Shaped  Tetrahedron  for  Each . 101 

Figure  23.  CDF  for  IDEAS  Generated  Mesh . 102 

Figure  24.  Relative  Error  in  Problem  Average  Flux  for  50%  Node  Variation . 104 

Figure  25.  Relative  Error  in  Problem  Average  Flux  for  50%  Node  Variation . 104 

Figure  26.  L,  Norm  of  Exiting  Face  Flux  for  50%  Variation . 105 

Figure  27.  LI  Norm  of  Exiting  Face  Flux  for  50%  Variation . 105 

Figure  28.  Average  Flux  for  Test  Problem  2 . Ill 

Figure  29.  Relative  Error  in  Average  Flux  for  Test  Problem  2 . Ill 

Figure  30.  Exiting  Face  Current  for  Test  Problem  2 . 112 

Figure  3 1 .  Relative  Error  in  Exiting  Face  Current  for  Test  Problem  2 . 113 

Figure  32.  Cut  Away  View  of  the  Test  Problem  3 

2511  Cell  Mesh . 117 

Figure  33 .  Problem  Average  Flux  for  Test  Problem  3 .  119 

Figure  34.  Cube  Material  Average  Flux  for  Test  Problem  3 . 120 

Figure  35.  Cylinder  Material  Average  Flux  for  Test  Problem  3 . 120 

viii 


Figure  36.  Exiting  Current  at  z  =  -1  cm  for  Test  Problem  3 . 121 

Figure  37.  Exiting  Current  at  y  =  -1  cm  for  Test  Problem  3 . 122 

Figure  38.  Problem  Average  Flux  for  Test  Problem  3 . 123 

Figure  39.  Exiting  Current  at  z  =  -1  cm  for  Test  Problem  3 . 124 

Figure  40.  Problem  Average  Flux  for  Test  Problem  4 . 126 

Figure  41.  Cube  Material  Average  Flux  for  Test  Problem  4 . 126 

Figure  42.  Cylinder  Material  Average  Flux  for  Test  Problem  4 . 127 

Figure  43.  Exiting  Face  Current  for  Test  Problem  4 . 128 

Figure  44.  Test  Problem  4  Exiting  Face  Flux  Distribution  for  the  EC  Method . 130 

Figure  45.  Test  Problem  4  Exiting  Face  Flux  Distribution  for  LC  Method . 131 

Figure  46.  Test  Problem  4  Exiting  Flux  Distribution  Computed  Using  EC 

on  Two  Coarsest  Meshes  Using  Nodes  Values  (Left  Side)  and  Using 

Internal  Data  Points  (Right  Side) . 132 

Figure  47.  Maximum  Pointwise  Relative  Error  of  the  Exiting  Face  Flux  Distribution  for 
Test  Problem  4 . 133 

Figure  48.  Average  Pointwise  Relative  Error  of  the  Exiting  Face  Flux  Distribution  for 

Test  Problem  4 . 134 

Figure  49  Details  of  the  number  of  single  moment  functions  that  have  to  be 
evaluated  to  compute  a  three  argument  function  by  recursion  on 
number  of  arguments . 148 

Figure  50.  Number  of  single-argument  moment  functions  that  have  to  be  computed  to 
evaluate  a  four-argument  function  completely  by  recursion  on  number  of 
arguments . 152 

Figure  5 1 .  Entering  Face  Showing  Node  Coordinates  in  Both  Local  uvw  and  in  Triangle 
Notation . 161 


IX 


Figure  52.  Area  in  (py,  px7)  Space  Containing  Valid  Solutions  to  the  2D  Rootsolve 

System . 166 

Figure  53.  Phase  Space  Volume  Containing  Solutions  to  3D  RootSolving  System.  .178 

Figure  54.  Case  1  Tetrahedron  in  Global  xyz  Coordinate  System  (left  figure)  and  in  Local 
uvw  Coordinate  System  (right  figure) . 183 

Figure  55.  Case  4  Tetrahedron  and  Resulting  Sub-Tetrahedra . 186 

Figure  56.  Case  4  Tetrahedron  Showing  Point  Labeling  Convention 

(Configuration  1) . 187 

Figure  57.  Case  4,  Sub-Tetrahedron  1  in  relation  to  Global,  Rotated, 

and  Local  Sub-Tetrahedron  Coordinate  Systems . 189 

Figure  58.  Case  4  Tetrahedron  in  Sub-Tetrahedron  1  Rotated  Coordinate  System.  .  .  192 

Figure  59.  Projection  of  Case  4,  Sub-Tetrahedron  1  Edges  into  the  v'w'-Plane . 192 

Figure  60.  Case  3  Tetrahedron  and  Resulting  Sub-Tetrahedra . 196 

Figure  61.  Case  3  Tetrahedra  and  Possible  Point  Labeling  Configurations . 197 

Figure  62.  Case -3  Tetrahedron  . 198 

Figure  63.  Side  View  of  Case  3  Tetrahedron . 201 

Figure  64.  Case  2  Tetrahedron  and  Resulting  Sub-Tetrahedra . 202 


x 


List  of  Tables 


Table  1.  Mesh  Sizes  Used  for  Test  Problem  1 . 83 

Table  2.  Execution  Time  per  Phase  Space  Cell  (seconds) . 106 

Table  3.  Comparison  of  One-Dimensional  Semi- Analytic  and  TETSN  Results  for 

Integral  Quantities  of  Test  Problem  2 .  115 

Table  4.  Volume  of  Mesh  Representation  of  Cylinder  (Source  Region)  for 

Test  Problem  3  and  Error  Relative  to  the  Exact  Volume  (0.78539  cnr) . 118 


xi 


AFIT/DS/ENP/96-0 1 


Abstract 

Characteristic  spatial  quadratures  for  discrete  ordinates  calculations  on  meshes  of 
arbitrary  tetrahedra  are  derived  and  tested,  including  the  step  (SC),  linear  (LC),  and 
exponential  (EC)  characteristic  quadratures  and  variants  that  assume  constant  distributions 
on  cell  faces.  Tetrahedral  meshes  accurately  model  curved  surfaces  with  few  cells.  A  split 
cell  approach  subdivides  tetrahedra  along  the  streaming  direction,  reducing  the  transport 
to  one  dimension.  Assumed  forms  of  the  cell  source  and  entering  flux  distributions  have 
sufficient  parameters  to  match  the  zeroth  and  first  spatial  moments.  These  parameters  are 
determined  by  analytically  inverting  a  linear  system  (LC),  or  by  numerical  inversion  using 
Newton's  method  (EC).  Efficient  algorithms  for  the  two-  and  three-dimensional 
rootsolves  are  derived.  The  constant  face  methods  proved  unacceptable  in  empirical 
testing.  Both  LC  and  EC  exhibited  third  order  convergence.  LC  provided  accurate 
results  on  cells  with  optical  thickness  on  the  order  of  one  mean  free  path  while  EC  was 
accurate  with  fewer,  thicker  cells.  LC  can  produce  negative  fluxes;  EC  is  strictly  positive. 
Although  more  costly  per  cell,  EC  is  robust  and  can  be  more  efficient  than  LC  or  SC  by 
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using  coarse  meshes. 


CHARACTERISTIC  SPATIAL  QUADRATURES 


FOR  DISCRETE  ORDINATES 
NEUTRAL  PARTICLE  TRANSPORT 
ON  ARBITRARY  TETRAHEDRAL  MESHES 

I.  Introduction 

Existing  codes  for  performing  three-dimensional  discrete  ordinates  calculations  are 
generally  limited  to  rectangular  parallelpiped  cells  and  low-order  spatial  quadratures.  The 
spatial  quadratures  currently  employed  in  these  codes  are  simple  to  implement  and  well 
understood,  but  are  also  very  inaccurate  unless  a  highly  refined  spatial  mesh  is  employed. 
Rectangular  parallelpiped  cells  are  especially  limiting  due  to  their  inability  to  adequately 
resolve  curved  boundaries  or  inclined  surfaces  without  excessive  refinement.  In  certain 
situations,  the  level  of  mesh  refinement  required  to  adequately  resolve  boundaries  may  be 
greater  than  that  required  to  ensure  accuracy  of  the  transport  solution  on  the  individual 
cells.  Non-rectangular  cells  defined  in  terms  of  curvilinear  coordinates  can  be  used  to 
alleviate  some  of  the  difficulties,  but  many  of  the  higher  order  solution  techniques  are  not 
amenable  to  derivation  in  curvilinear  coordinates.  This  research  addresses  these 
limitations  by  deriving  and  implementing  highly  accurate  characteristic  spatial  quadratures 
for  meshes  of  arbitrary  tetrahedra. 


1.1.  Background 

The  fundamental  equation  describing  the  neutral  particle  transport  process  is  the 
linear  Boltzmann  transport  equation  (BTE)  (Lewis  and  Miller,  1993:  1) .  The  BTE 
describes  the  balance  of  neutral  particles  in  a  seven-dimensional  phase  space:  three  spatial 
dimensions,  two  angular  dimensions,  energy,  and  time  (Lewis  and  Miller,  1993:  24-42). 
The  BTE  relates  the  time  rate  of  change  of  the  neutral  particle  density  at  a  point  r,  moving 
with  energy  E,  in  direction  Q,  at  time  t,  to  the  difference  between  sources  and  losses: 


—  -J-+d'V  +  o(F ,E,t) 
v  at 


ty(r,E,Q,t)=S(r,E,Q,t)  , 


(1) 


where  the  angular  flux  density,  <]r,  is  the  flux  of  particles  at  point  r,  with  energy  E,  moving 
in  direction  Q,  at  time  t,  the  bracketed  term  accounts  for  particle  losses  via  streaming  and 
collisions,  and  the  right  hand  side  accounts  for  all  sources  of  particles.  The  equation  is 
examined  in  detail  in  section  2.1.  Analytic  solutions  to  the  BTE  exist  only  for  highly 
idealized  cases.  In  general,  solutions  are  obtained  using  approximate,  numerical 
techniques.  In  certain  cases,  the  diffusion  equation,  which  can  be  derived  from  the  BTE 
using  several  simplifying  assumptions  about  the  angular  dependence,  can  be  used  to  obtain 
approximate  solutions.  Methods  for  solving  the  fully  angular  dependent  BTE  fall  into  two 
broad  categories:  Monte  Carlo  and  deterministic. 
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1.1.1  Diffusion  Approximation 

Accounting  for  the  angular  dependence  of  the  transported  particles  introduces 
considerable  complexity  and  represents  a  large  part  of  the  computational  effort  in  most 
solution  methods.  For  many  applications  this  added  complexity  is  unnecessary  and 
reasonable  results  can  be  achieved  by  assuming  that  the  flux  is  only  weakly  dependent  on 
angle,  i.e.  that  the  flux  is  linearly  anisotropic  (Duderstadt,  1976: 133-140).  Under  this 
assumption,  the  BTE  can  be  used  as  a  starting  point  in  the  derivation  of  the  diffusion 
equation: 


-  V-Z)(^£)V<f>(7,£,0  +  o  (r,E)$(r,E,t) 

v  dt 

(2) 

=  fdE' ,t)  +  S(r,E,t)  . 


Where  D,  the  diffusion  coefficient,  determines  the  rate  at  which  the  diffusion  occurs,  and 
the  scalar  flux,  rf>,  is  defined  as 

§{r,E,t)  =  pQi|;(rA£,r)  .  (3) 

Since  the  assumption  of  weak  angular  dependence  works  well  in  the  cores  of 
thermal  nuclear  reactors,  the  diffusion  approximation  is  widely  used  in  the  nuclear 
industry.  However,  the  application  of  diffusion  theory  requires  that  the  migration  process 
be  dominated  by  scatter,  and  that  the  process  be  far  removed  from  any  sources  or  sinks, 
i.e.  large  gradients  which  introduce  an  angular  dependence.  However,  for  some  very 
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important  applications,  such  as  near  control  rods  in  thermal  reactors,  blankets  in  fast 
reactors,  or  deep  penetration  shielding  type  problems,  the  diffusion  approximation  is 
inadequate,  and  transport  theory  must  be  used  (Duderstadt,  1979: 138). 

1.1.2.  Monte  Carlo  Methods 

Monte  Carlo  is  a  technique  by  which  individual  particle's  paths  are  simulated. 

These  histories  are  accumulated  and  statistics  are  formed  to  estimate  desired  processes. 
One  advantage  of  Monte  Carlo  is  its  ability  to  handle  complex  three-dimensional 
geometries.  Another  advantage  of  Monte  Carlo  is  the  continuous  treatment  of  space  and 
angle  variables  which  eliminates  the  discretization  errors  associated  with  deterministic 
methods.  Monte  Carlo  errors  take  the  form  of  stochastic  uncertainties  (Lewis  and  Miller, 
1993:  296).  One  disadvantage  of  the  technique  is  that  the  accuracy  of  the  estimators  is 
dependent  on  statistical  variance.  Reducing  the  variance  may  require  a  large  number  of 
samples,  which  is  computationally  expensive.  Some  variance  reduction  techniques  are 
available  and  improve  the  accuracy  of  estimators  while  reducing  the  computational 
expense.  However,  use  of  these  techniques  is  as  much  art  as  science  (Lewis  and  Miller, 
1993:  297-298). 

1.1.3.  Deterministic  Methods 

Unlike  Monte  Carlo,  which  can  be  understood  without  reference  to  the  BTE, 
deterministic  methods  solve  the  BTE  in  one  form  or  another  by  using  various  techniques 
to  treat  the  angular  and  spatial  integrations.  The  most  widely  used  method  for  solving  the 
BTE  is  discrete  ordinates  (Lewis  and  Miller,  1993:  1 16).  In  discrete  ordinates,  the  angular 
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variable  is  discretized,  and  the  BTE  is  solved  for  a  discrete  set  of  directions,  or  ordinates. 
The  angular  integrations  then  reduce  to  weighted  sums.  The  spatial  variable  is  also 
discretized,  and  treated  using  an  appropriate  spatial  quadrature.  The  discrete  ordinates 
method  was  first  utilized  to  perform  transport  calculations  in  stellar  atmospheres 
(Chandrasekhar:  1960).  The  method  was  adapted  for  use  in  the  neutron  transport 
community  by  Carlson  and  others  working  at  Los  Alamos  National  Laboratory  (Carlson, 
2386:  1958;  Carlson  and  Lathrop:  1968) 

The  basic  computational  unit  in  a  discrete  ordinates  calculation  is  the  spatial  cell, 
and  as  such,  computational  cost  and  storage  requirements  are  directly  proportional  to  the 
number  of  cells  in  the  calculation.  Initial  discrete  ordinates  codes  utilized  the  diamond 
difference  (DD)  spatial  quadrature  which  uses  a  finite  difference  approximation  to  the 
derivative  term  in  the  BTE  along  with  several  auxiliary  assumptions  (Lewis  and  Miller, 
1993:  128-131).  The  diamond  difference  method  is  easy  to  implement  and 
computationally  inexpensive,  but  can  be  very  inaccurate  on  cells  with  optical  thickness 
greater  than  about  a  tenth  of  a  mean  free  path  (mfp).  On  thicker  cells  (in  slab  geometry) 
and  in  two  or  three  spatial  dimensions,  the  method  may  generate  unphysical  negative 
fluxes.  In  two-dimensional  rectangular  geometry,  unphysical  oscillations  in  the  DD 
generated  solutions  have  been  observed  (Petrovic  and  Haghighat:  1995).  The 
shortcomings  of  diamond  difference  prompted  efforts  to  develop  more  robust  schemes 
leading  to  the  development  of  the  linear  discontinuous  (Lewis  and  Miller,  1993:  134-135), 
exponential  (Barbucci  and  Di  Pasquantonio:  1977),  linear  nodal  (Walters  and  O'Dell, 
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1981)  and  linear  characteristic  (Larsen  and  Alcouffe,  1981)  methods.  The  linear  nodal  and 
linear  characteristic  methods  provide  accurate  solutions  on  cells  with  optical  thickness  on 
the  order  of  1  mfp  (generally  2-4  mfp),  but  can  still  produce  negative  fluxes  on  thicker 
cells. 

The  desire  to  have  good  coarse  mesh  accuracy  and  strict  positivity  prompted  the 
development  of  the  nonlinear  characteristic  methods:  step  adaptive  (SA),  linear  adaptive 
(LA)  and  most  recently,  exponential  characteristic  (EC).  Each  of  the  methods  solve  the 
BTE  on  a  spatial  cell  by  integrating  along  a  characteristic.  Assumed  forms  of  the  source 
distribution  are  used  having  parameters  that  are  chosen  to  match  the  spatial  moments  of 
the  distribution  on  the  spatial  cell.  The  SA,  LA,  and  EC  methods  use  increasingly 
smoother  source  distribution  approximations.  The  SA  method  assumes  a  piecewise 
constant  source  distribution  (Mathews:  1991;  Mathews  and  Minor:  1993),  and  the  LA 
method  assumes  a  piecewise  linear  cell  source  distribution  (Mathews  and  Minor:  1993). 
The  exponential  characteristic  method  (Sjonden:  1992;  Mathews,  Sjonden  and  Minor: 
1994;  Walters  and  Wareing:  1994;  Minor:  1993;  Mathews  and  Brennan:  1995)  uses  an 
exponential  source  distribution  which  can  be  thought  of  as  an  extension  of  the  S  A  and  LA 
distributions  to  a  C00  function.  All  of  the  nonlinear  methods  have  the  desirable  property  of 
being  strictly  positive.  The  EC  method,  in  addition  to  being  strictly  positive  exhibits 
excellent  accuracy  on  optically  thick  cells  (optical  thickness  >  10  mfp). 

Due  to  the  relatively  high  computational  cost  per  spatial  cell  of  the  EC  method, 
efforts  have  begun  to  develop  related  methods  that  use  exponential  source  and  flux 
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representations,  but  are  less  computationally  intensive.  These  methods  include  the 
exponential  discontinuous  (ED)  method  (Walters  and  Wareing:  1996)  and  the  nonlinear 
comer  balance  method  (NLCB)  (Castrianni  and  Adams:  1995).  Both  of  these  methods  are 
relatively  new  and  as  of  the  present  time,  their  performance  relative  to  the  EC  method  has 
not  been  firmly  established. 

Central  to  the  development  of  improved  solution  methods  has  been  the  application 
of  the  methods  to  unstructured  grids.  Monte  Carlo  methods  are  readily  adapted  to  general 
three-dimensional  problems,  making  them  the  method  of  choice  for  three-dimensional 
transport  calculations.  While  there  are  many  situations  for  which  Monte  Carlo  is  the 
appropriate  choice,  there  are  problems  for  which  discrete  ordinates  could  be  more 
efficient,  particularly  when  spatial  distributions  are  required.  The  current  generation  of 
three-dimensional  discrete  ordinates  codes  use  diamond  difference  or  linear  nodal  spatial 
differencing  (Alcouffe:  1995;  Rhoades  and  Azmy:  1995)  and  meshes  of  rectangular 
parallelpipeds.  Rectangular  parallelpiped  meshes  are  especially  limiting  for  several 
reasons.  First,  many  small  cells  are  required  to  adequately  resolve  curved  boundaries  or 
inclined  surfaces.  Further,  if  the  mesh  needs  to  be  refined  to  provide  higher  resolution  in  a 
small  region  of  the  problem,  it  must  be  refined  everywhere.  Decreasing  the  cell  size  by  a 
factor  of  n  in  each  direction  results  in  an  m  increase  in  the  total  number  of  cells,  and  a 
corresponding  nJ  increase  in  computation  cost  and  storage. 

General  geometry  methods  for  performing  efficient  neutral  particle  transport  in 
multiple  dimensions  are  beginning  to  be  developed.  There  have  been  several  recent  efforts 
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to  develop  methods  for  solving  the  BTE  on  unstructured  meshes  of  cells  in  two- 
dimensions.  Comer-balance  schemes  (Adams:  1993;  Castrianni  and  Adams:  1995)  solve 
the  BTE  on  an  unstructured  grid  of  polygons.  The  step  and  linear  characteristic  methods 
have  been  extended  to  2D  triangular  (Miller:  1993)  and  arbitrary  polygonal  (Groves  and 
Pevey:  1995)  meshes,  and  the  exponential  characteristic  method  has  also  been  applied  to 
unstructured  meshes  of  triangles  (Mathews  and  Brennan:  1995). 

The  goal  of  this  research  is  to  derive  and  implement  the  characteristic  methods  on 
meshes  of  arbitrary  (unstructured)  tetrahedra.  This  will  address  the  need  for  a  general 
geometry,  three-dimensional  discrete  ordinates  transport  code  that  may  provide  a  more 
computationally  efficient  alternative  to  Monte  Carlo  and  the  existing  three-dimensional 
discrete  ordinates  codes. 

The  characteristic  methods  allow  for  accurate  solutions  on  large,  or  optically  thick, 
cells  (cell  through  which  the  streaming  length  is  many  mean  free  paths)  by  integrating  the 
transport  equation  along  streaming,  or  characteristic,  lines  of  particle  flow  (characteristic 
lines  of  the  BTE).  One  drawback  of  the  method  is  the  fact  that  the  characteristic 
integrations  must  be  analytically  evaluated  to  derive  the  methods.  In  curvilinear 
coordinates,  the  characteristics  are  not  straight  lines,  and  the  resulting  integrations  are 
intractable  in  those  coordinate  systems.  However,  curvilinear  coordinates  are  useful  only 
for  problems  having  the  required  symmetries.  For  example,  cylindrical  coordinates  work 
well  for  a  cylinder,  or  a  set  of  concentric  cylinders,  but  are  useless  for  two  cylinders  with 
distinct  axes. 
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Using  meshes  of  arbitrarily  oriented  tetrahedra  addresses  many  of  the  previously 
mentioned  difficulties  associated  with  discrete  ordinates  solution  methods.  First,  they  are 
defined  in  terms  of  a  Cartesian  coordinate  system,  so  the  characteristics  are  straight  lines, 
and  the  resulting  integrations  can  be  analytically  evaluated.  Second,  curved  surfaces  that 
required  the  use  of  a  large  number  of  parallelpiped  cells  to  model  can  be  reasonably 
approximated  using  a  much  smaller  collection  of  arbitrary  tetrahedra.  Finally,  global  mesh 
refinements  are  not  necessary.  Mesh  refinements  in  localized  regions  of  the  problem  can 
be  accomplished  with  minimum  impact  on  the  size  of  the  cells  in  the  rest  of  the  mesh. 

1.2.  Purpose 

The  purpose  of  this  research  is  to  derive,  implement,  and  evaluate  the  step,  linear, 
and  exponential  characteristic  methods  for  meshes  of  arbitrary  tetrahedra.  In  addition, 
variants  of  the  linear  and  exponential  characteristic  methods  using  constant  boundary  flux 
representations  will  also  be  examined.  Dr.  W.  F.  (Wally)  Walters  suggested  investigating 
the  constant  face  methods  to  determine  if  the  higher  order  source  distribution 
approximation  would  increase  the  efficiency  of  the  step  characteristic  method  (Walters, 
1994). 

Several  of  the  required  tools  were  brought  forward  from  previous  efforts.  The 
algorithm  used  to  test  the  methods  is  a  three-dimensional  adaptation  of  Miller's  (1993) 
arbitrary  triangle  algorithm  for  the  LC  method  on  meshes  of  arbitrary  triangles.  In 
deriving  the  EC  method  for  2D  rectangular  meshes.  Minor  (1993)  developed  routines  to 
numerically  compute  exponential  moment  functions  of  all  orders  with  up  to  three 


9 


arguments.  For  the  tetrahedral  mesh  methods,  routines  to  compute  four-argument 
moment  functions  were  developed,  and  all  of  the  multi-dimensional  functions  routines 
were  rewritten  to  increase  their  efficiency. 

Specifically,  the  implementation  of  the  quadratures  requires: 

a.  Extending  the  triangle  discrete  ordinates  algorithm  (Miller,  1993)  to  three 
dimensions. 

b.  Developing  definitions  of,  appropriate  coordinate  systems  for,  and  optimal 
means  of  passing  spatial  moments  on  tetrahedral  faces. 

c.  Developing  methods  for  splitting  arbitrary  tetrahedra  along  the  particle 
streaming  direction  and  generating  the  sub-tetrahedra  and  necessary 
transformations  between  the  parent  and  sub-tetrahedra  coordinate  systems. 

d.  deriving  equations  for  computing  the  cell  and  exiting  flux  moments  for  each  of 
the  methods  in  the  appropriate  coordinate  systems. 

e.  Developing  efficient  algorithms  for  determining  coefficients  of  the  entering 
flux  and  source  distributions  using  knowledge  of  the  zeroth  and  first  spatial 
moments. 

f.  Developing  efficient  routines  to  evaluate  multi-dimensional  exponential 
moment  functions. 

g.  Developing  a  method  to  generate  the  required  tetrahedral  meshes  and  create 
appropriate  data  structures  for  input  to  the  algorithm. 

1.3.  Scope 

Although  with  straightforward  modification,  the  methods  derived  here  can  be 
applied  to  a  wide  range  of  problems,  the  derivations  include  only  enough  complexity  to 
demonstrate  the  effectiveness  of  the  spatial  quadratures  on  meshes  of  arbitrary  tetrahedra. 
First,  only  the  time-independent,  or  steady  state  case  is  considered.  Time  dependent 
problems  are  beyond  the  scope  of  this  research. 
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As  the  spatial  quadratures  essentially  invert  the  streaming-collision  operator  on 
the  left  hand  side  of  the  BTE,  the  complete  details  of  the  source  term  on  the  right  hand 
side  of  the  equation  need  not  be  considered  to  demonstrate  the  methods.  The  derivations 
assume  that  all  scatter  is  isotropic.  Although  the  addition  of  anisotropic  scatter  to  the 
method  is  straightforward  (Lewis  and  Miller,  1993:61-77),  its  inclusion  is  not  necessary  to 
demonstrate  correct  operation  of  the  methods  and  was  therefore  omitted.  Similarly,  only 
non-multiplying,  or  no  fission,  cases  were  considered.  Sub-critical  multiplication  simply 
requires  the  inclusion  of  an  additional  term  on  the  right  hand  side  of  the  equation,  and 
critical  and  super-critical  systems  require  both  a  slightly  altered  form  of  the  BTE  and  the 
modified  source  term,  but  neither  requires  any  modification  to  the  spatial  quadratures 
(Lewis  and  Miller,  1993:90-103). 

No  type  of  acceleration  technique  (Lewis  and  Miller,  1993:  145-150)  was 
examined  for  accelerating  the  convergence  of  the  methods  for  optically  thick  problems 
with  large  scattering  to  total  cross  section  ratios.  There  are  many  problems  for  which  the 
methods  as  derived  are  entirely  applicable.  At  this  point  it  is  not  clear  if  the  any  of  the 
methods  need  to  be  modified  in  any  way  to  accommodate  an  accelerator,  but  the  ability  to 
handle  problems  with  cross  sections  approaching  the  diffusion  limit  is  beyond  the  scope  of 
the  current  effort.  It  is  anticipated,  however,  that  a  linearization  acceleration  technique 
would  be  applicable  to  the  EC  method  developed  here  ( Wareing,  Walters,  and  Morel, 
1995). 
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1.4.  Organization 

The  remainder  of  this  document  consists  of  five  chapters  and  four  appendices. 
Chapter  II  discusses  the  discrete  ordinates  method,  identifies  the  role  of  the  spatial 
quadrature,  and  introduces  the  general  form  of  the  characteristic  spatial  quadratures. 
Chapter  III  develops  the  characteristic  methods  for  arbitrary  tetrahedral  meshes,  including 
development  of  coordinate  systems,  the  splitting  of  arbitrary  tetrahedra  to  produce  special 
case  tetrahedra,  the  special  case  moment  equations  for  each  of  the  methods  derived,  and 
the  derivation  and  uses  of  the  balance  equations.  Chapter  IV  presents  the  solution 
algorithm,  the  details  of  the  iterative  solution  procedure,  and  concludes  with  a  pseudocode 
outline  of  the  entire  code.  Chapter  V  presents  the  results  of  numerical  tests  of  the  code  on 
several  test  problems  and  reports  on  the  basic  operating  characteristics  of  the  code 
including  convergence  rates,  sensitivity  to  mesh  variations,  speed,  and  accuracy.  The  final 
chapter  presents  my  conclusions  and  recommendations  for  further  efforts. 

There  are  five  appendices  that  provide  a  more  detailed  discussion  of  procedures 
used  in  the  algorithm.  Appendix  A  details  the  numerical  algorithm  used  to  evaluate  the 
exponential  moments  functions.  Appendix  B  deals  with  the  derivation  and  operation  of 
the  two-  and  three-dimensional  root-solving  routines.  Appendix  C  provides  further  details 
of  the  tetrahedral  splitting  routines.  Appendix  D  contains  a  discussion  of  the  use  of  the 
balance  equations  to  compute  cell  flux  moments.  Appendix  E  discusses  mesh  generation 
using  the  IDEAS  (Lawry,  1991)  code. 
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II.  Derivation  of  the  Characteristic  Spatial  Quadratures 

2.1.  The  Boltzmann  Transport  Equation 

The  BTE,  introduced  in  Section  1.1,  describes  the  balance  of  neutral  particles  in  a 
seven-dimensional  phase  space.  The  BTE  relates  the  time  rate  of  change  of  the  neutral 
particle  density  at  a  point  r,  moving  with  energy  E,  in  direction  Q,  at  time  t,  to  the 
difference  between  sources  and  losses: 

- i|r (?, E, Cl ,  t)  +  [□ • ' V  +  o (r,E,  t) ]  i|r  (r ,E, Q ,  t)  =  S(r ,E,Cl ,  t) ,  (4) 

v  at 

where  the  angular  flux  density,  t|r,  is  the  flux  of  particles  at  point  r,  with  energy  E,  moving 
in  direction  Q,  at  time  t.  The  bracketed  term  on  the  left-hand  side  of  the  equation  is 
termed  the  streaming-collision  operator.  The  streaming-collision  operator  consists  of  the 
particle-streaming-loss  and  the  collision-loss  terms.  The  total  interaction  cross  section,  a, 
represents  the  probability  at  point  r,  per  unit  path  length  of  particle  travel,  that  a  particle 
having  energy  E  will  undergo  a  collision.  The  total  cross  section  term  in  the  streaming 
collision  operator  accounts  for  losses  caused  by  both  scatter  and  absorption,  o  =  os+oa. 
Particles  with  a  given  energy  and  direction,  removed  by  scatter,  become  a  source  of 
particles  in  other  energies  and  directions.  Particles  lost  by  absorption  can  also  result  in 
sources  to  other  energies  and  directions  via  secondary  emission  from  fission,  (n,  2n), 

(n,  y),  (n,  a)  and  other  similar  reactions. 

The  right  hand  side  of  the  equation  contains  the  sum  of  all  sources  of  particles,  and 
can  be  expanded  into 
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S(r,E,€l,t)  =  qs{r,E,Q,t)+qf(r,E,Cl,t)+qext(r,E,Cl,t)  , 


(5) 


which  contains  contributions  from  scatter,  qs,  fission,  qf,  and  other  (extrinsic)  sources, 

fiext- 

The  scattering  source  term,  qs,  accounts  for  the  production  of  particles  with 
energy  E  moving  in  direction  Q  as  a  result  of  scattering  collisions. 


qs(r,E,Q,t)  = 


/ 


dQ'oi^E'-E^Cl'-Qjy^i^E'^j)  . 


(6) 


Equation  (6)  states  that  the  source  of  particles  with  energy  E,  moving  in  direction  Q,  at 
point  r,  and  time  t,  is  simply  the  sum  of  all  particles  scattered  into  energy  E  and  direction 
Q,  from  other  energies  E'  and  directions  Q'. 

The  fission  term,  qf,  accounts  for  the  production  of  particles  (specifically 
neutrons)  with  a  given  energy  and  direction  due  to  fission, 

qf{r,E,Cl,t)  =  x(E)fdE' favour, E',t)  t |/(f,£/,Q/,0  .  (7) 

The  term  x(E)  accounts  for  the  fact  that  the  energy  distribution  of  fission  neutrons  is 
independent  of  the  energy  of  the  particle  that  caused  the  fission.  Since  number  of 
neutrons  produced  is  the  desired  quantity,  it  is  conventional  to  use  vcy,  the  product  of  the 
fission  cross  section,  of,  and  the  average  number  of  neutron  emitted  per  fission,  v.  Note 
that  the  fission  cross  section  contains  no  reference  to  the  incident  or  exiting  neutron 
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directions.  This  is  because  fission  reactions  produce  neutrons  with  an  isotropic  angular 
distribution,  and  the  probability  of  fission  is  independent  of  the  direction  of  motion  of  the 
incident  neutron. 

The  extrinsic  source  term,  qext,  accounts  for  particles  produced  by  means  other 
than  scattering  or  fission.  This  may  include  material  sources,  such  as  particles  produced  as 
a  result  of  spontaneous  fission  or  decay,  cosmic  rays,  etc. 

Although  the  BTE  is  solved  for  tjr,  the  angular  flux  itself  is  not  a  physically  useful 
quantity.  The  desired  products  are  usually  particle  currents  or  reaction  rates.  A  more 
useful  quantity  for  computing  reaction  rates  is  the  scalar  flux,  cj>,  defined  by 

<p{f,E,/)  =  JdQ\}r{f,E,Q, t)  .  (8) 

Whereas  the  angular  flux  distribution,  ijr  dfdE  dCl ,  can  be  thought  of  as  the  total  path 
length  per  unit  time  traveled  by  particles  in  the  phase  space  volume  {dfdE  dCl )  moving 
in  a  given  direction,  Q,  the  scalar  flux  distribution,  <j)  dfdE,  is  just  the  total  path  length 
rate  of  particles  in  the  phase  space  volume  traveling  {dfdE)  in  all  directions.  Because 
reaction  cross  sections  are  almost  always  independent  of  the  incident  particle  direction, 
reaction  rates  are  the  product  of  the  reaction  cross  section  and  the  scalar  flux. 

Particle  leakage  rates  are  governed  by  currents.  Again,  one  is  generally  interested 
in  the  total  number  of  particles  penetrating  a  shield  or  leaking  through  a  crack,  not  in  the 
specific  angular  distribution  of  the  penetrating  particles.  The  current  vector  is  defined  as 
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(9) 


■/ 


J(r,E,t)=  dQ  Q  \|r(f, E,Q,t) 


The  net  number  of  particles  per  unit  time  per  unit  energy  per  unit  surface  area  crossing  a 
surface  traveling  in  a  direction  with  a  component  in  the  fi  (the  outward  surface  normal) 
direction  is  then  just 

Jn(r,E,t)  =  j'dQthQ  i| s(r,E,Q,t)  =  n- J(r,E,t)  .  (10) 


As  is  obvious  from  the  equations,  computing  the  angular  flux  correctly  is  essential  to 
computing  the  scalar  flux  and  current  correctly. 

2.2.  Discrete  Ordinates  Formulation 

The  time-independent  form  of  the  BTE  in  Cartesian  coordinates  is 


+  A  -^-+^+a(x,y,z,E) 
ox  dy  dz 


ty{x,y,z,E,Cl)  =  S(x,y,z,E,Q ) 


01) 


where  p  =  ex-Q  ,  r|  =ey-Q  ,  and  £  =  e/Q  are  the  direction  cosines  of  Q.  Defining 
the  differential  solid  angle,  dQ  as 


_  d(j>  J0sin0  _  dw  dE, 
In  2  2 n~  ’ 


(12) 


gives  the  normalization. 
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(13) 


(Lewis  and  Miller,  1993,  28).  Assuming  a  non-multiplying  system  and  isotropic  scatter, 
and  using  equations  (8)  and  (13),  the  source  term  reduces  to 


S(x,y,z,E)  =  JdE/as(x,y,z,E/^E)<t>(x,y,z,E/)  +qea(<x,ytz,E) 


(14) 


Rather  than  solving  over  a  continuous  angular  variable,  the  discrete  ordinates  method 
solves  the  BTE  for  a  discrete  set  of  directions,  or  ordinates.  The  discretized  (in  angle) 
form  of  Equation  (1 1)  for  a  single  direction,  is 


u  —  +  ri  —  +  £  —  +  o(x,y,z,E ) 
indy  "dz  K  ' 


ty„(x,y,z,E)  =  S(x,y,z,E)  . 


(15) 


The  angular  integration  used  to  compute  the  scalar  flux  and  current  is  now  approximated 
by  a  weighted  sum  of  the  form 


$(x,y,z,E )  =  wn  t| r„(x,y,z,E)  .  (16) 

n  -  l 

The  set  of  direction  cosines,  (pn,  r|n,  ^n)  and  weights,  wn,  is  termed  the  angular 
quadrature  set.  Quadrature  sets  are  chosen  both  to  assure  accuracy  of  the  numerical 
integrations  and  to  provide  required  symmetries  for  reflective  boundary  conditions.  The 
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set  of  angular  quadratures  that  are  invariant  under  90  degree  rotations  about  any  axis  are 
called  the  level  symmetric  quadratures.  If  the  symmetry  requirement  can  be  relaxed,  other 
quadratures,  tailored  to  the  specifics  of  the  problem,  can  be  derived  (Lewis  and  Miller, 
1993:  162-165;  Abu-Shumays,  1977). 

The  directions  associated  with  a  given  quadrature  set  also  determine  the  directions 
in  which  the  particles  are  allowed  to  travel.  Allowing  particle  flow  in  only  a  discrete  set  of 
directions  can  lead  to  unphysical  behaviors  of  the  flux,  called  ray  effects  (Lewis  and 
Miller,  1993:  194-203),  as  described  in  section  2.7. 

The  energy  dependence  is  also  handled  by  discretization.  In  this  case,  the  energy 
spectrum  is  divided  into  a  set  of  discrete  bins,  or  energy  groups.  The  discretization 
produces  a  set  of  equations  (15)  for  each  energy  group.  These  equations  are  coupled 
through  the  scattering  term  because  as  particles  lose  (or  gain)  energy  through  scatter,  they 
are  transported  from  group  to  group.  Discretizing  equation  (1 5)  in  energy  yields 


(17) 


The  group  fluxes,  t|;^g  and  4>n  g,  are  given  by 


7 

A in,g(x,y,z)  =  J  dE'^n(x,y,z,E’)  , 


(18) 
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where  the  limits  of  integration  are  the  lower  and  upper  bounds  of  energy  group  g.  The 
group  source  term  is  given  by 

g-i 

Sg(x,y,z)  =  E  4>g'(x,y,z)  (2Q) 

+  o^g(x,y,z)  <| \{x,y,z)  +  q^g(x,y,z)  , 

where  the  group  scattering  cross  section,  agUg,  is  the  probability  that  a  particle  in  energy 
group  g1  undergoing  a  collision  will  emerge  in  energy  group  g.  The  group  subscripts 
increase  for  decreasing  energy:  group  1  is  the  highest  energy  group;  group  G  is  the  lowest. 
Equation  (20)  only  allows  for  down-scatter,  i.e.  particles  may  not  gain  energy  in  collisions. 
For  many  problems,  including  photon  transport  and  fast  reactors,  this  assumption  is 
strictly  valid  (Lewis  and  Miller,  1993,  79).  For  thermal  reactors,  there  is  only  appreciable 
up  scatter  for  neutrons  with  thermal  energies.  In  this  case  equation  (20)  can  be  used  if 
the  energy  range  in  which  up  scatter  occurs  is  contained  in  a  single  energy  group 
(Duderstadt,  1976,  291-295). 

Spatial  discretization  is  accomplished  by  dividing  the  domain  into  a  set  of  spatial 
cells.  The  cross  sections  are  assumed  constant  within  each  spatial  cell.  There  is  then  an 
equation  (17)  for  each  spatial  cell  for  each  angle  and  energy  group.  Once  discretized,  the 


solution  is  obtained  using  a  spatial  quadrature  to  invert  the  streaming-collision  operator  in 
equation  (17). 

2.2.1  Boundary  Conditions 

To  solve  the  BTE,  the  flux  distribution  entering  the  solution  volume,  V,  across  the 
bounding  surface,  T,  must  be  known.  Specifically,  if  n  has  an  outward  orientation  with 
respect  to  T,  then  xj/n  (r)  for  n  •  Qn  <  0  and  re  T  are  required  as  boundary  data.  For 
a  known  incident  flux,  a(r),  the  boundary  condition  is 

=  Qn- n<0 ,  reT  (21) 

A  vacuum  boundary  condition  results  when  T  (r)  =  0 . 

The  surface  source  and  vacuum  boundary  conditions  are  explicit.  Implicit 
boundary  conditions,  can  be  used  to  take  advantage  of  symmetries.  An  albedo  boundary 
condition  is  given  by 

ynAr)  =  aoty,o{r),  Q-w<0, fel\  (22) 

where  aa  is  the  albedo  and  Q  is  the  angle  resulting  from  specular  reflection  of  particles 
traveling  in  the  Qn  /  direction.  The  special  case  where  a  =  1  is  termed  a  symmetry 
boundary.  Other  types  of  implicit  boundaries  include  the  white  boundary,  in  which  all 
particles  are  reflected  back  from  a  surface  with  an  isotropic  (Lambertian)  distribution,  and 
the  periodic  boundary  condition,  where  the  distribution  entering  through  one  boundary  is 
set  equal  to  the  distribution  leaving  through  the  corresponding  boundary  in  a  periodic 
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lattice  (Lewis  and  Miller,  1993:  25-26).  The  current  implementation  of  the  code  supports 
both  types  of  explicit  boundary  conditions  and  the  albedo  boundary  condition. 

2.3.  Discrete  Ordinates  Solution 

The  coupled  set  of  equations  (17),  (19),  and  (20),  together  with  the  appropriate 
boundary  conditions,  are  solved  using  an  iterative  procedure  termed  the  iteration  on  the 
scattering  source.  The  procedure  consists  of  an  outer  iteration  through  each  of  the  energy 
groups,  and  in  inner  iteration  for  each  energy  group  that  consists  of  a  loop  through  each  of 
the  angles  in  the  quadrature,  and  a  loop  through  each  of  the  cells  in  the  spatial  mesh.  The 
outer  iteration  is  the  iteration  on  the  between  group  scattering  source  in  equation  (19). 

The  inner  iteration  is  the  iteration  on  the  within  group  scattering  source. 

The  inner  iteration  consists  of  a  “walk”  through  the  spatial  mesh  for  each  angle 
until  the  within  group  scattering  source  for  the  current  energy  group  is  converged.  The 
angle  determines  the  direction  of  particle  flow,  and  hence  the  direction  of  walk  through 
the  spatial  mesh.  Equation  (17)  is  solved  for  each  cell  in  the  spatial  mesh  using  a  spatial 
quadrature,  defined  in  detail  in  Section  2.4.  Input  information  to  the  spatial  quadrature 
consists  of  spatial  moments  of  the  fluxes  on  the  entering  faces  of  the  cell  (determined  by 
the  current  streaming  direction)  and  the  source  within  the  cell.  Initially,  the  cell  source  is 
determined  using  an  initial  value  for  the  group  scalar  fluxes  in  equation  (19)  (usually  zero) 
to  compute  the  group  source  for  each  cell  in  the  spatial  mesh.  Given  the  source  and 
entering  flux  information,  the  cell  interior  and  exiting  fluxes  are  computed  using  the  spatial 
quadrature.  The  exiting  fluxes  for  the  current  cell  then  become  the  entering  fluxes  for  the 


21 


neighboring  downstream  cells.  The  procedure  continues  in  the  direction  of  particle  flow 
until  each  of  the  cells  in  the  spatial  mesh  have  been  computed.  The  current  angle 
contribution  to  the  scalar  flux  is  now  known  and  added  to  a  running  total  using  equation 
(16).  A  new  angle  is  then  chosen  and  the  procedure  repeated  until  all  the  angles  in  the 
quadrature  set  have  been  used.  At  this  point,  each  term  in  the  sum  of  equation  (16)  has 
been  accumulated.  The  current  group  scalar  fluxes  for  each  cell  are  compared  to  those 
from  the  previous  iteration.  If  the  relative  change  in  values  is  larger  than  some 
predetermined  tolerance,  the  source  term  for  each  cell  is  updated  using  the  new  group 
scalar  flux  in  equation  (19),  the  current  iteration  scalar  flux  value  is  saved  and  used  for 
comparison  at  the  end  of  the  next  iteration,  and  the  iteration  procedure  repeated.  If  the 
relative  change  in  values  is  less  than  some  predetermined  tolerance,  the  method  is  said  to 
have  converged.  The  current  group  scalar  flux  is  then  used  in  equation  (19)  to  update  the 
down-scatter  source  for  the  lower-energy  groups.  This  completes  an  inner  iteration.  This 
procedure  is  then  repeated  for  each  of  the  energy  groups,  completing  the  outer  iteration. 
An  outline  of  the  discrete  ordinates  iterative  procedure  follows: 

Determine  initial  values  for  the  group  scalar  fluxes 
For  each  energy  group 
Do: 

Update  scattering  source  with  previous  iteration  scalar  fluxes 
For  each  angle: 

For  each  spatial  cell: 
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Compute  cell  and  exiting  angular  flux  data 


Pass  cell  exiting  flux  data  to  downstream  cells  as 
input  information 

Update  scalar  flux  sum  with  current  cell  angular 
flux 

Next  Spatial  Cell 
Next  angle 

Check  for  convergence  by  comparing  previous  iteration  scalar  flux 
to  current  value 

Loop  until  converged 

Update  lower  energy  group  scattering  sources  with  current  group  scalar 
flux 

Next  Energy  Group 

Convergence  of  the  method  can  be  explained  by  examining  the  physical 
interpretation  of  the  iteration  on  the  scattering  source.  Initially,  the  scattering  source  is 
zero  (or  some  initial  guess),  and  all  particles  in  the  problem  are  either  present  initially, 
arise  from  internal  sources,  or  enter  through  the  problem  boundaries.  Since  none  of  these 
particles  have  scattered,  the  flux  computed  in  the  first  iteration  is  termed  the  first-flight,  or 
uncollided,  flux.  For  the  second  iteration,  the  source  consists  of  particles  arising  from 
extrinsic  sources  together  with  a  contribution  from  particles  that  have  been  scattered  once. 
The  third  iteration  source  contains  contributions  from  particles  that  have  been  scattered  at 
most  twice,  etc.  If  the  physics  of  the  problem  is  such  that  the  average  particle  undergoes 
only  a  few  scatters  before  either  being  absorbed  or  leaking  out  of  one  of  the  boundaries. 
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then  successive  contributions  to  the  scattering  source  from  particles  that  have  undergone 
one  more  scatter  should  ultimately  decrease  to  the  point  where  the  updated  angular  fluxes 
are  not  significantly  different  from  the  previous  iteration  fluxes.  At  this  point  the  method 
is  said  to  have  converged. 

2.4.  Role  of  the  Spatial  Quadrature 

The  spatial  quadrature  numerically  inverts  the  streaming-collision  operator  of 
equation  (17)  on  a  computational  cell.  The  inputs  to  the  spatial  quadrature  are  the  spatial 
moments  of  the  scattering  source  distribution  within  the  cell  and  the  flux  distribution  on 
the  entering  face(s)  of  the  cell.  The  spatial  quadratures  derived  here  use  assumed  analytic 
forms  of  the  input  distributions  with  parameters  chosen  to  match  either  the  zeroth  or  both 
the  zeroth  and  first  moments  of  the  distributions.  Using  these  parameters,  the  quadrature 
computes  the  moments  of  the  flux  distribution  within  the  cell  volume  and  on  the  exiting 
face(s). 

The  spatial  moments  are  defined  as  weighted  integrals  over  either  the  cell  face  or 
cell  volume.  The  zeroth  moment  is  simply  the  average,  and  the  weighting  function  is  one. 
The  first  moment  describes  how  the  flux  or  source  is  distributed  with  respect  to  one  of  the 
coordinate  axes,  and  is  weighted  accordingly.  Moments  of  the  entering  flux  distribution 
are  defined  as 


(23) 
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where  ijrm  (s,,  s2)  is  the  assumed  form  of  the  entering  flux  distribution,  and 

r(?i  ,s2)=ea+sle1+s2e2,  (24) 

where  e0  is  a  convenient  point  on  the  face,  such  as  a  vertex,  and  ex  and  e2  are 
orthogonal  unit  vectors  that  span  the  plane  of  the  face.  The  moments  of  the  source 
distribution  are  defined  as 

S(x  pXj,x3),  (25) 

where  S(xl5  x2,  x3)  is  the  assumed  form  of  the  source  distribution,  and  xl5  x2,  and  x3 
parameterize  the  cell  volume. 

The  actual  values  of  the  entering  flux  moments  are  known,  either  from  the 
boundaiy  conditions  or  passed  from  a  upstream  cell.  The  values  of  the  source  distribution 
moments  are  also  known,  either  from  an  initial  guess  or  computed  using  previous  iteration 
results.  Equations  (24)  and  (25)  are  evaluated  analytically  giving  systems  of  equations  in 
these  known  moments  having  the  parameters  of  the  assumed  distributions  as  unknowns. 
The  parameters  are  then  determined  by  inverting  the  system,  either  by  direct  analytical  or 
numerical  inversion  for  linear  systems,  or  by  a  multi-dimensional  root  solve  for  non-linear 
systems. 
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The  spatial  quadrature  inverts  equation  (17),  giving  an  expression  for  the  cell  flux 
distribution  in  terms  of  the  source  and  entering  flux  parameters.  This  procedure  is 
described  in  detail  in  the  following  section.  Once  an  expression  for  the  flux  in  the  cell  is 
known,  it  can  be  used  to  determine  the  moments  of  the  cell  and  exiting  flux  distributions. 
The  moments  of  the  cell  flux  distribution  are  defined  as 
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and  the  moments  of  the  exiting  flux  distribution  are  defined  as 


(  1  out') 

A  A 

,  out 

Vsl 

■iJK 

.  out 

+*  J 

A*) 

t  [*i  (st ,  s2)  ,x2  ( Jj  ,s2),x3  (5j  ,s2)  ]  (27) 


where  i|;ceI1(x1,  x2,  x3)  is  the  cell  flux  distribution  given  in  terms  of  the  entering  flux  and 
scattering  source  parameters,  and  Sj  and  s2  are  defined  for  the  exiting  face.  The  cell 
angular  flux  moments  are  used  in  equation  (19)  to  generate  the  scalar  flux  moments  which 
are  then  used  to  update  the  scattering  source  distribution  for  the  next  iteration: 
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where  G  is  the  total  number  of  energy  groups,  the  subscript  g  indicates  the  group  number, 
and  the  superscript  k  indicates  the  iteration  number  (Note  that  equation  (28)  assumes 
down-scatter  only).  The  exiting  flux  moments  become  the  entering  flux  moments  for 
downstream  cells. 

Using  the  details  of  the  spatial  quadrature,  the  discrete  ordinates  iterative 
algorithm  is 

Initialize  scalar  fluxes 
For  each  energy  group: 

Do: 

Update  source  moments  using  previous  iteration  scalar  fluxes 

Determine  the  scattering  source  distribution  parameters  for  each 

spatial  cell  using  known  source  moments 

For  each  angle: 

For  each  spatial  cell: 

Determine  the  parameters  of  the  entering face  flux 
distributions  from  the  known  moments 

Compute  cell  and  exiting  flux  moments  using  the 
entering  and  source  distribution  parameters 
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Pass  cell  exiting  flux  moments  to  downstream  cells 
as  input  information 

Accumulate  current  angle  contribution  to  scalar 
flux  sum 

Next  spatial  cell 
Next  Angle 

Check  for  within  group  convergence  by  checking  relative  error 
between  current  and  previous  iteration  scalar  fluxes 

Loop  Until  Converged 

Update  lower  energy  group  scattering  sources  using  current  group  scalar 
fluxes 

Next  energy  group 

2.5.  Moments  of  the  Boltzmann  Transport  Equation 

The  moments  integrals  can  also  be  applied  to  the  discretized  BTE  to  produce 
conservation  or  balance  equations.  The  conservation  equations  describe  the  balance  of 
particles  over  the  spatial  cell.  Existing  spatial  quadratures  utilize  either  only  the  zeroth  or 
the  zeroth  and  first  moments.  Methods  that  use  only  the  zeroth  moment  satisfy  the  zeroth 
moment  balance  equation,  methods  that  use  both  zeroth  and  first  moments  satisfy  both  the 
zeroth  and  first  moment  balance  equation.  In  general,  methods  that  utilize  first  moments 
exhibit  higher  orders  of  convergence  and  are  more  accurate  than  methods  that  use  only  the 
zeroth  moment.  However,  the  cost  of  increased  accuracy  and  faster  convergence  is 
increased  computational  effort  and  storage  cost  for  each  cell  (Lewis  and  Miller,  1993: 
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135).  Nevertheless,  the  more  complicated  methods  can  be  more  efficient  by  producing  the 
same  accuracy  with  fewer  cells. 

2.6.  Characteristic  Solutions  to  the  Boltzmann  Transport  Equation 

The  angular  and  energy  dependence  of  the  BTE  are  handled  by  iterating  on  a 
discretized  form  of  the  equation.  The  spatial  dependence  is  treated  by  use  of  the  spatial 
quadrature.  Characteristic  spatial  quadratures  solve  the  discretized  form  of  the  BTE, 
Equation  (17),  by  integrating  along  a  characteristic,  or  particle  streaming  direction.  The 
general  form  of  the  characteristic  solution  for  a  single  ordinate  has  the  form 


tJrCnQJ  =  iJr(F0,QJe 


J'ds/e~a’(s~s')S(s/), 

0 


(29) 


where  s  is  the  position  along  the  characteristic  line  starting  at  r0  traveling  in  the  Qn 
direction,  and  r  =  r0  +sQn  (Alcouffe,  R.E.  and  E.W.  Larsen,  1981:3-6).  The  different 
characteristic  methods  are  distinguished  by  the  assumed  form  of  the  entering  flux 
distribution,  i|/(r0,Qn),  and  source  distribution,  S(s). 

2.6.1.  Step  Characteristic 

The  step  characteristic  (SC)  method  assumes  that  the  entering  flux  and  source 
distributions  are  constant. 


S(x,y,z)=SA 

x r(Sl,s2)  =  ^, 


(30) 


29 


where  Sj  and  Sj  are  parameters  defining  points  on  the  entering  face.  The  SC  method  is 
computationally  inexpensive  compared  to  higher  order  methods.  The  method  satisfies 
only  the  zeroth  moment  balance  equation.  The  method  has  been  demonstrated  to  be  first- 
order  convergent  in  slab  geometry  (Larsen  and  Miller,  1980:  78)  and  in  two-dimensional 
rectangular  and  triangular  geometries  (Minor,  1993:  89-95;  Miller,  Mathews  and  Brennan, 
1996). 

The  representation  of  the  entering  flux  and  source  distributions  as  constants  can 
lead  to  inaccuracies  caused  by  numerical  diffusion.  The  inaccuracies  and  low  order  of 
convergence  resulting  from  the  assumption  of  constant  entering  flux  and  source  suggest 
the  use  of  a  higher  order  approximation  to  the  distributions. 

2.6.2.  Linear  Characteristic 

The  linear  characteristic  (LC)  method  models  the  source  and  entering  flux  as  linear 
distributions  within  the  cell  and  on  the  entering  face: 


S(x,y,z)  =  a  +  bx  +  cy  +  dz , 
ilr'"(51,52)  =  oc+P151  +p252. 


(31) 


The  LC  method  has  been  extensively  evaluated  and  has  been  shown  to  perform  well  on  a 
wide  range  of  problems  (Larsen  and  Miller,  1980:  106-107;  Mathews,  Sjonden,  and 
Minor,  1994:33-36;  Miller,  1993:  70-113;  Minor,  1993:  83-139) .  The  method  conserves 
both  zeroth  and  first  moments,  and  has  been  shown  to  be  fourth  order  convergent  in  slab 
geometry  (Larsen  and  Miller,  1980:  80)  and  third  order  convergent  in  two  dimensional 
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rectangular  and  triangular  geometries.  (Larsen  and  Alcouffe,  1981;  Miller,  Mathews  and 
Brennan,  1996). 

The  linear  representation  of  the  source  in  the  spatial  cell  can  lead  to  problems  for 
optically  thick  cells,  i.e.  cells  that  are  more  than  a  few  mean-ffee-paths  thick.  If  the  source 
distribution  is  such  that  the  linear  approximation  is  negative  on  the  exiting  face  of  the  cell, 
an  unphysical  negative  flux  can  result  (Minor,  1993:9-12).  Should  this  occur,  it  can  be 
remedied  in  one  of  three  ways.  First,  the  cell  size  can  be  reduced,  thereby  reducing  the 
optical  thickness  of  all  cells.  This  approach  increases  the  number  of  cells  required  to 
model  the  geometry,  and  increases  the  computational  cost.  A  second  solution  is  to  invoke 
a  negative  flux  fix  up.  There  are  several  procedures  available  for  fixing  up  negative  fluxes, 
but  each  destroys  first  moment  conservation,  degrading  the  accuracy  of  the  solution 
(Larsen  and  Alcouffe,  1981:  101-106).  At  third  approach  is  to  abandon  the  linear 
approximation  and  choose  a  distribution  that  is  strictly  positive  and  more  realistically 
models  the  actual  distribution. 

2.6.3.  Exponential  Characteristic 

The  exponential  characteristic  (EC)  method  models  both  the  entering  flux  and 
source  distributions  as  exponentials, 

S(x,y,z)=aebx+cy+dz 

&s+Qs  (32) 
and  i|f”(Sj,s2)  =  ae  '  1  2  . 
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In  optically  thick  cells  without  embedded  sources,  the  exponential  distribution  provides  a 
reasonable  approximation  to  the  actual  distributions.  Further,  the  exponential  distribution 
has  the  added  advantage  of  being  strictly  positive.  The  method  conserves  both  zeroth  and 
first  moments,  and  displays  fourth  order  convergence  in  slab  geometry  (Mathews, 

Sjonden,  and  Minor,  1994:  33-37)  and  third  order  in  convergence  in  two  dimensional 
rectangular  (Minor,  1993:  89-95)  and  triangular  (Mathews  and  Brennan,  1995)  geometry. 
The  method  has  been  demonstrated  to  perform  exceptionally  well  in  deep  shielding  type 
problems,  where  the  distributions  are  strongly  attenuated  across  the  domain  of  the 
problem  (Mathews,  Sjonden,  and  Minor,  1994;  Walters  and  Wareing,  1994;  Minor  and 
Mathews,  1995;  and  Mathews  and  Brennan,  1995).  When  initially  introduced,  the  method 
was  much  less  computationally  efficient  per  spatial  cell  than  existing  quadratures 
(approximately  20  times  the  cost  of  LC),  but  its  accuracy  on  coarse  meshes  allowed  it  to 
obtain  highly  accurate  results  using  a  relatively  small  number  of  cells  resulting  in  less 
overall  computational  effort  and  greatly  reduced  storage  requirements  (Mathews,  Sjonden, 
and  Minor,  1993;  Minor,  1993;  Mathews  and  Brennan,  1996).  However,  recent 
improvements  in  the  implementation  and  special  function  algorithms  have  reduced  the  cost 
to  about  twice  that  of  the  LC  method  per  spatial  cell  on  triangular  meshes  (Mathews  and 
Brennan,  1996). 

2.6.4.  LC  and  EC  Variants  -  Constant  Face  Methods 

Variants  of  the  first  moment  conserving  methods  (LC,  EC)  can  be  derived  in  which 
the  source  distribution  is  expanded  in  a  linear  or  exponential  function  but  the  face  flux  is 
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assumed  constant.  These  methods  are  also  derived,  implemented,  and  tested  as  part  of 
this  project. 

2.7.  Discretization  Errors 

Discretization  of  the  space  and  angle  variables  invariably  introduces  errors.  In  two 
or  more  dimensions,  systematic  errors  caused  by  angular  discretization  can  occur.  In 
problems  with  small  scattering-to-total  cross  section  ratios,  the  neutron  distribution  may 
contain  peaks  that  coincide  with  the  streaming  directions  of  the  particles  that  are 
inconsistent  with  the  physical  realities  of  the  problem.  These  non-physical  oscillations  are 
known  as  ray  effects.  Ray  effects  can  be  reduced  by  several  means  including  increasing 
the  order  of  the  quadrature  or  reverting  to  a  product  quadrature  (Lewis  and  Miller,  1993: 
195-203;  Abu-Shumays,  1977). 

Spatial  discretization  can  also  contribute  to  a  non-physical  redistribution  of 
particles  known  as  numerical  diffusion.  Numerical  diffusion  is  caused  by  smearing  the 
input  flux  distribution  over  the  entire  face  of  the  cell  or  of  the  source  distribution  over  the 
interior  of  the  cell,  either  of  which  could  result  in  an  erroneous  redistribution  of  particles 
throughout  the  problem.  In  some  cases,  the  numerical  diffusion  and  ray  effects  can,  to  an 
extent,  mask  each  other.  Particles  incorrectly  bunched  along  a  streaming  direction  due  to 
ray  effects  can  to  some  extent  be  fortuitously  redistributed  almost  correctly  by  the  effects 
of  the  numerical  diffusion,  nevertheless,  the  solution  accuracy  is  degraded  by  both  effects. 
Numerical  diffusion  can  be  mitigated  by  carrying  first  moments  of  the  distribution  along 
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the  cell  faces.  This  helps  to  maintain  the  proper  distribution  of  the  angular  flux  as  it  is 
passed  between  cells  (Mathews,  1983:  VI-5  -  VI-8). 


III.  Derivation  of  the  Characteristic  Methods 


for  Arbitrary  Tetrahedral  Meshes 

3.1.  Application  of  the  Characteristic  Methods  to  Arbitrary  Tetrahedral  Meshes 

Trying  to  solve  the  transport  equation  on  a  mesh  of  arbitrary  tetrahedra,  where 
everything  is  defined  in  terms  of  a  single,  global  coordinate  system  introduces  prohibitive 
complexity.  In  essence,  under  such  a  scheme  every  integration  used  to  compute  cell 
moments  must  be  considered  as  a  special  case.  It  is  possible  to  choose  a  local  coordinate 
system  for  each  tetrahedron  that  standardizes  the  form  of  the  source  and  flux  distribution 
moments.  This,  however,  complicates  the  characteristic  solution  because  the 
characteristics  are  arbitrary  in  the  local  coordinate  system,  and  the  characteristic 
integration  has  to  be  treated  as  a  special  case  for  each  tetrahedron.  Miller  (1993)  solved 
these  difficulties  in  two  dimensions  by  developing  a  procedure  whereby  arbitrary  triangles 
were  split  along  the  characteristic  direction.  The  tetrahedral  mesh  algorithm  presented 
here  is  a  direct  extension  of  Miller's  technique  to  three  dimensions. 

Consider  a  tetrahedron  oriented  in  space  in  such  a  way  that  the  particle  streaming 
direction  is  coincident  with  one  of  the  edges,  as  shown  in  Figure  1 .  The  local  u,v,w 
coordinate  system  is  oriented  such  that  the  u-axis  is  coincident  with  the  streaming 
direction  and  the  base  of  the  tetrahedron  lies  in  the  uv  plane.  We  term  this  special  case 
tetrahedron  a  case  1  tetrahedron,  referring  to  the  fact  that  there  is  only  one  inflow  face  and 
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one  outflow  face.  For  a  case  1  tetrahedron  in  the  given  coordinate  system,  the  one-group 
form  of  equation  (17)  reduces  to  the  one-dimensional  form. 


d_ 

du 


+  o 


\lr(u,v,w)=S(u,v,w) , 


(33) 


where  =  eu,  \in  =  Qn-eu  =  1 ,  x\n  =  Q„-ev  =  0  and ,  5  =  =  0  and  we  assume  that 

the  cross  sections  are  constant  within  the  cell.  Because  the  characteristic  direction  is 
parallel  to  the  u-axis,  the  characteristic  equation,  equation  (29),  becomes 


t|f(M,v,w)  =  t|/;.n[w;n(v,w),v,w]c 


du/e-a(u-u')S(u/,v,w)  .  (34) 


Substitution  of  the  assumed  form  of  the  entering  flux  and  source  distribution  into  equation 
(34)  gives  a  functional  form  for  the  cell  flux  distribution  in  terms  of  the  entering  flux  and 
source  distribution  coefficients.  This  can  then  be  used  to  compute  the  flux  distribution 
moments  in  the  cell  volume  and  on  the  cell  outflow  face. 

3.2.  Splitting  of  Arbitrary  Tetrahedra  into  Sub-Tetrahedra 

The  algorithm  developed  here  uses  a  splitting  routine  that  divides  arbitrarily- 
oriented  tetrahedra  along  the  streaming  direction  such  that  a  number  of  sub-tetrahedra  are 
produced,  each  of  which  has  exactly  one  input  and  one  output  face.  Each  such  subcell  is  a 
case  1  tetrahedron  because  the  edge  at  the  intersection  of  the  remaining  two  faces  is 
oriented  along  the  streaming  direction.  Using  this  method,  a  local  coordinate  system  can 
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be  erected  for  each  sub-tetrahedron  identical  to  that  shown  in  Figure  1 .  The  characteristic 
equation  is  then  identical  for  all  sub-tetrahedra  and  only  a  single,  consistent  set  of  analytic 
equations  for  the  angular  flux  moments  are  required;  no  special  cases  need  be  considered. 
This  simplifies  both  the  derivation  of  the  analytic  form  of  the  moments  and  the 
computation  of  the  numerical  values  of  the  moments. 

The  possible  results  of  tetrahedral  splitting  are  shown  in  Figure  2.  Each  possibility 
is  labeled  by  the  resulting  number  of  case  1  tetrahedra  resulting  from  the  split.  In  reality, 
there  are  two  general  cases,  case  4  and  case  3,  and  two  degenerate  cases,  case  2  and  case 
1 .  If  the  streaming  direction  is  incident  on  the  tetrahedron  such  that  there  are  two  inflow 
and  two  outflow  faces,  the  tetrahedron  is  case  4.  The  case  4  tetrahedron  is  split  along  a 
line  parallel  to  the  streaming  direction  that  enters  one  edge  and  exits  the  opposite  edge 
resulting  in  four  case  1  tetrahedra. 

When  three  sub-tetrahedra  result,  the  tetrahedron  is  either  case  3  or  case  -3.  A 
case  3  tetrahedron  is  characterized  by  having  three  input  faces  and  one  output  face.  A 
case  -3  tetrahedron  has  one  input  face  and  three  output  faces.  A  case  3  tetrahedron  is  split 
along  a  line  parallel  to  the  streaming  direction  that  enters  through  an  apex  and  exits 
through  the  opposite  (outflow)  face.  A  case  -3  tetrahedron  is  split  along  a  line  parallel  to 
the  streaming  direction  that  enters  through  a  face  and  exits  through  the  apex  opposite  the 
inflow  face. 

A  case  2  tetrahedron  has  two  entering  faces  and  one  exiting  face.  Likewise,  a  case 
-2  tetrahedron  has  one  entering  face  and  two  exiting  faces.  In  either  case  there  is  one  face 
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parallel  to  the  streaming  direction.  A  case  2  tetrahedron  is  split  along  a  line  parallel  to  the 
streaming  direction  that  enters  through  an  apex  and  exits  the  face  opposite  the  apex.  A 
case  -2  tetrahedron  is  split  along  a  line  that  enters  through  a  face  and  exits  through  the 
opposite  apex.  Either  case  can  be  seen  to  be  a  degenerate  case  3  tetrahedron  with  the 
entering  or  exiting  point  on  the  face  located  on  an  edge. 

Appendix  C  contains  a  more  detailed  discussion  of  each  case  and  provides  the 
details  needed  to  implement  the  splitting  procedure. 


Figure  1 .  Case  1  Tetrahedron  in  the  Local  uvw  Coordinate  System, 
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Case  -2 


Case  2 


Case  1  Tetrhaedra 
one  input  face 
one  output  face 


Case  2  Tetrahedron 
two  input  faces 
one  output  face 
Case  -2  Tetrahedron 
one  input  face 
two  output  faces 


Case  3  Tetrahedron 
three  input  faces 
one  output  face 
Case  -3  Tetrahedron 
one  input  face 
three  output  faces 


Case  4  Tetrahedron 
two  input  faces 
two  output  faces 


Figure  2.  Tetrahedra  Splitting  and  Resulting  Number  of  Case  1  Sub-Tetrahedra. 
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3.3.  Definition  of  “Arbitrary”  Tetrahedra 

When  referring  to  meshes  of  arbitrary  tetrahedra,  the  term  “arbitrary”  refers  to  the 
sizes,  shape,  and  orientation  of  the  tetrahedron  in  space  and  the  characteristics  of  the 
individual  tetrahedra.  There  is  no  restriction  on  how  individual  tetrahedra  are  oriented 
with  respect  to  the  defining  coordinate  system.  There  are  no  size  or  aspect  ratio  limits  for 
the  individual  tetrahedra.  Any  number  of  tetrahedra  may  share  an  individual  node.  The 
only  restriction  is  that  each  face  must  be  shared  by  only  two  tetrahedra,  i.e.  faces  may  not 
be  split  so  that  a  tetrahedron  has  more  than  a  single  neighbor  across  any  given  face. 

3.4.  Coordinate  Systems 

Like  its  triangle  predecessor,  the  arbitrary  tetrahedra  algorithm  computes  such 
quantities  as  source  distribution  coefficients,  entering  flux  distribution  coefficients,  and  cell 
flux  moments  in  local  (orthogonal)  coordinate  systems  designed  to  simplify  both  the 
derivation  and  resulting  form  of  the  moment  equations.  There  are  five  separate  coordinate 
systems  used  in  the  algorithm:  global  (x,y,z),  centroid  (xc,yc,zc),  face  (xF,yF),  source 
(u',v',w'),  and  local  sub-tetrahedron  (u,v,w). 

The  global  coordinate  system  is  generated  with  the  problem  geometry.  All 
tetrahedra  nodes  are  specified  in  the  global  xyz  coordinate  system.  The  global  system 
remains  fixed  and  independent  throughout  the  problem  execution.  Each  (parent)  cell  in 
the  mesh  has  its  own  local,  centroid  coordinate  system.  For  a  given  tetrahedron,  the 
centroid  system  is  generated  by  translating  the  origin  of  the  global  xyz 
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system  to  the  tetrahedron  centroid.  All  cell  averaged  quantities  are  computed  in  the 
centroid  system.  The  relation  between  the  global  and  centroid  systems  is  shown  in 
Figure  3. 


The  face  coordinate  system,  shown  in  Figure  4,  is  defined  for  each  face  in  the  mesh 
without  regard  to  the  tetrahedron  with  which  the  face  is  associated;  i.e.  the  common  face 
between  two  adjacent  tetrahedra  has  a  single  coordinate  system  defined  over  it.  This 
trivializes  the  passing  of  moments  between  adjacent  tetrahedra:  they  are  identical.  For  a 
given  face,  the  system  is  defined  so  that  the  xF  axis  lies  along  an  edge  and  the  apex  not  on 
the  xF  axis  has  a  positive  yF  coordinate.  The  zF  axis  is  normal  to  the  face,  in  the  direction 
necessary  to  make  a  right-handed  coordinate  system.  The  selection  of  the  apex  defining 
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the  origin  and  the  edge  defining  the  xF  axis  are  arbitrary,  however  once  chosen,  they 
remain  fixed  throughout  the  problem  execution.  Moments  of  the  entering  and  exiting 
angular  flux  distribution  are  computed  in  terms  of  the  face  system  for  each  tetrahedron. 
Note  that  if  the  xF  axis  is  chosen  to  be  the  longest  edge  of  the  triangle,  0  <  y  <1,  which 
may  provide  better  numerical  conditioning  in  the  computation  of  the  face  distribution 
coefficients. 

The  source  coordinate  system  is  a  local  system  defined  for  each  original 
tetrahedron  in  the  mesh,  and  is  shown  in  Figure  5.  It  is  used  in  the  first  order  methods  as  a 


convenient  system  for  computing  the  coefficients  of  the  source  distribution.  The 
orientation  of  the  tetrahedron  in  the  source  system  is  identical  to  the  orientation  of  the 
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Case  1  tetrahedron  in  the  local  sub-tetrahedron  coordinate  system  shown  in  Figure  1.  The 
choice  of  the  apex  defining  the  system  origin  and  the  edge  defining  the  u-axis  are  arbitrary, 
but  once  chosen,  they  remain  fixed  throughout  the  problem  execution.  Note  that  if  the  u- 
axis  is  chosen  to  run  along  the  longest  edge  of  the  tetrahedron,  then  Uj  >  u2  >  0  and  Uj  > 
u3  >  0.  This  may  provide  better  numerical  conditioning  in  the  computation  of  the  source 
coefficients. 


Figure  5.  Coordinate  System  Used  to  Compute  the  Coefficients  of  the 
Source  Distribution. 


The  local  sub-tetrahedron  uvw  coordinate  system  was  introduced  in  the  previous 
section,  and  is  shown  in  Figure  1 .  The  spatial  quadrature  is  applied  in  the  local  coordinate 
system,  i.e.  the  source  and  entering  flux  distribution  coefficients  are  transformed  into  the 
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local  uvw  system,  and  the  cell  and  exiting  flux  moments  are  computed  with  respect  to  the 
local  system  for  each  case  1  tetrahedron. 

3.5.  Mathematical  Preliminaries:  The  Exponential  Moments  Functions 

The  exponential  moment  functions  are  a  class  of  functions  originally  developed  to 
deal  with  removable  singularities  encountered  in  discrete  ordinates  methods.  Numerical 
problems  arise  when  the  equations  experience  catastrophic  cancellation  as  they  approach 
their  removable  singularity.  For  example 

limjt'-^Ul,  (35) 

but  when  evaluated  numerically  as  x-0,  all  precision  is  lost  due  to  catastrophic 
cancellation  in  the  numerator.  When  these  types  of  singularities  are  encountered,  the 
expression  must  be  evaluated  as  a  special  case  to  preserve  accuracy. 

The  first  generation  of  moment  functions  (Walters,  1986:  192-194;  Walters,  1986: 
22-24),  were  moment  functions  of  a  single  argument.  Mathews  (Mathews,  Sjonden,  and 
Minor,  1995:  27)  later  developed  a  general  integral  definition  of  the  functions  and  showed 
that  it  implied  the  recursion  on  order  relationship.  Minor  (Minor,  1993:  46-59)  extended 
the  definition  to  include  multiple  arguments,  and  termed  this  broader  class  the  exponential 
moment  functions.  In  addition.  Minor  developed  numerical  routines  to  evaluate  moment 
functions  with  up  to  three  arguments.  All  but  the  single  argument  routines  were  rewritten 
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as  part  of  this  project  to  increase  their  computational  efficiency.  In  addition,  routines  were 
developed  to  compute  moment  functions  with  four  arguments. 

Minor's  dissertation  (Minor,  1993)  contains  a  complete  description  of  the  functions 
and  the  single  argument  algorithms.  To  facilitate  both  the  derivation  of  the  methods  and 
the  explanation  of  the  moments  function  algorithm,  this  section  begins  with  an  overview  of 
the  functions. 

3.5.1.  Definition  and  Relationships 

The  exponential  moment  functions  are  defined  as 

1  “l  Um- 1 

Mn(xx,x2,...,xJ=  ^dux^du2...  j  dum(l  -i/)n  e~x‘Ul  e(x'~X2)U2 ...  e(Xm'rXm)Um ,  (36) 

0  0  0 

where  n  is  termed  the  order  of  the  function,  and  the  dimensionality  is  determined  by  the 
number  of  arguments,  m.  Using  this  definition,  an  equation  for  the  zeroth  order  function 
of  a  single  argument  is 

M0(x)  =  i — .  (37) 


Given  the  above,  higher  order  functions  of  a  single  argument  can  be  found  by  recursion  on 
order  using 


Mn(x)  = 


1 


X 


(38) 
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Recursion  relationships  also  exist  for  functions  of  more  than  one  variable,  for  example, 
any  function  of  more  than  one  variable  can  be  constructed  of  lower  dimension  functions 
using  recursion  on  dimension,  given  by 


Mn  (Xj  ,x2, ...  ,xm  Mn  (X2 ,  x3 , . . . ,  xm ) 


x  -  X, 


(39) 


Functions  of  more  than  one  argument  with  order  greater  than  0  can  also  be  generated 
using  recursion  on  order  given  by 


)  = 


3^2  ’ ,Xm  - 1  )  i  (Xj  ,X2, ...  ,Xm) 


(40) 


Equation  (39)  is  used  in  the  numerical  routines  to  compute  moments  functions  when  the 
values  of  the  arguments  are  far  enough  apart  to  avoid  numerical  difficulties.  Equation 
(40)  is  used  in  the  moments  functions  routine  to  compute  higher  order  functions  having 
the  same  arguments. 

The  ratio  of  the  first  to  zeroth  order  moment  functions  turns  out  to  be  a  useful 
quantity  for  computing  the  flux  moments  in  the  EC  method,  hence,  for  simplicity  of 
notation  we  define 


P(*i,*2> 


My  (Xj  ,X2,  ...  ,Xm) 

M0(xi,x2,...,xm) 


(41) 
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3.6.  Derivation  of  the  Characteristic  Methods  for  Meshes  of  Arbitrary  Tetrahedra 

This  section  examines  the  derivation  and  application  of  the  characteristic  spatial 
quadratures  for  arbitrary  tetrahedra.  The  spatial  walk  consists  of  the  application  of  the 
spatial  quadrature  to  each  of  the  tetrahedra  in  the  mesh.  As  the  steps  involved  in  applying 
the  quadrature  are  exactly  the  same  for  each  cell,  we  need  only  examine  one.  We  begin 
with  a  discussion  of  the  relatively  simplistic  SC  method,  explaining  the  steps  required, 
derivation  of  the  equations,  and  coordinate  system  transformations.  The  discussion  then 
proceeds  through  the  LC  and  EC  methods,  each  of  which  introduces  further  complexity. 
After  introducing  the  standard  versions,  we  examine  simplified  variations  of  the  LC  and 
EC  methods  in  which  the  source  distributions  retain  their  original  form  but  the  entering 
flux  distribution  is  assumed  constant. 

3.6.1.  Step  Characteristic  Method  (SC) 

The  first  step  in  developing  a  characteristic  scheme  is  to  specify  the  functional 
form  of  the  source  and  entering  flux  distributions.  For  the  SC  method,  both  are  assumed 
constant  and  equal  to  the  average  value  of  the  distribution.  The  source  distribution 
average  is  simply  the  zeroth  spatial  moment  computed  using  equation  (24).  The  entering 
flux  average  moment  is  given  either  by  the  boundary  conditions  or  passed  from  the 
adjacent  upstream  cell(s).  The  averages  are  invariant  with  respect  to  coordinate 
transformations,  hence  they  are  identical  for  each  sub-tetrahedron. 

Given  the  source  and  entering  flux  moments,  the  tetrahedron  is  split  along  the 
streaming  direction  producing  a  number  of  sub-tetrahedra  consistent  with  the  case  of  the 
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parent.  Local  coordinate  systems  are  then  established  for  each  subcell.  The  node 
coordinates  are  transformed  from  the  global  to  the  local  uvw  coordinate  system  using 


/  \ 
u 

(  \ 
X-Xx 

V 

=  R 

y  ~y  i 

z-z. 

\  V 

(42) 


where 


R  = 


/a  a  a  a  a  A  ^ 

e  -  e  e  ~e  e  •  *? 

u  x  u  y  u  z 


e  -e  e  -e  e  -e 

V  X  V  y  V  Z 

A,  /V  ^ 

.  e  -e  e  me  e  me  . 

\  w  x  w  y  w  z ) 


(43) 


and  Xj,  yl9  and  z}  are  the  coordinates  of  the  uvw  origin  in  the  global  coordinate  system. 

The  moments  of  the  subcell  flux  distribution  are  computed  in  the  local  coordinate 
system  by  substituting  the  constant  forms  of  the  source  and  entering  flux  distributions  into 
equation  (34)  to  get  a  functional  form  for  ij/(u,v,w)  used  in  equation  (26).  The  resulting 
integral  can  be  separated  into  the  source  and  entering  flux  contributions  and  the  average 
flux  moment  is  just  their  sum: 


■  subcell  _  ,  subcell  .  subcell 
VA  ”  VA ,S  +  ™Ajn 


(44) 


The  source  contribution  to  the  subcell  flux  moment  is  defined  as 
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subcell 

]A,S 


ds  SAe'^u-s)  / 


subcell 


udv’w'> 


(45) 


where  the  subcell- volume  integral  operator,  ( •  )subcell,  is  defined  as  (reference  Figure  1) 


vR(w)  u„Jv^w) 


0  VZ.(W)  udv’w'> 


where 


X  x  out  out 

Uout(V’W)  =  Cl  V+C2  W+Ul 

,  x  in  in 

W/«(V’W)  =  C  1  V+C7  W  * 


vL(w)  =  —  w  ,  and  , 


w. 


(v-  _v?) 

VR (W)  =  - W  +  V2  ■ 


(47) 


Additionally, 


out  U2  U\ 


out  _  V2Ul~V2Ul+UlV3~U2V3 


Cj"  =  —  ,  and  , 


in  W3V2"W2V3 


C2  ~ 


V2W3 


(48) 


Evaluating  equation  (45)  gives 


¥AUTU  =  ulSAM,(e) 


(49) 


where  e  =  o  ux  is  the  cell  optical  thickness. 

The  entering  flux  contribution  to  the  subcell  flux  moment  is  defined  as 


% 


subcell 
A,  in 


< 


subcell 


which,  when  evaluated  gives 


,  subcell  ~  .ju  w 

'I'm*  =  3  i| ^M2(€) 


(50) 


(51) 
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The  moments  of  the  exiting  flux  distribution  are  computed  in  the  same  way; 
substitute  the  constant  forms  of  the  source  and  entering  flux  distributions  into  equation 
(34)  to  get  i|/(u,v,w)  and  use  the  result  in  equation  (27).  The  total  exiting  flux  moment 
can  then  be  written  as  the  sum  of  the  source  and  exiting  flux  contributions: 


,  out  ,  out 

=%,S 


out 

A ,  in  ' 


(52) 


The  source  contribution  to  the  exiting  flux  moment  is  defined  as 


/  “ou,(V’W ) 

U=  \  J  dsSAe 


^ ^  face 


“*,(  V,w) 


(53) 


where  the  face  integral  operator,  ( •  )face,  is  defined  as 


(54) 


Performing  the  integration  yields  the  source  contribution: 


0«i^M2(e). 


(55) 


The  entering  flux  contribution  to  the  exiting  flux  moments  is  defined  as 
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(56) 


out 
'A,  in 


( 


_  -°[“0„,(v»w)-«Jv.w)] 
va  e 


which,  when  evaluated  gives  the  entering  flux  contribution: 

1|C  =  211^(6).  (57) 


When  the  subcell  moments  are  computed,  they  must  be  assembled  to  generate  the 
parent  cell  results.  As  already  stated,  the  average  moments  are  invariant  with  respect  to 
coordinate  transformations,  hence  they  have  the  same  value  in  the  centroid  and  global  face 
systems  as  in  the  uvw  system,  and  no  transformation  is  necessary.  The  parent  cell  average 
flux  is  simply  the  volume  weighted  average  of  the  subcell  average  fluxes: 


1  ^  subcells 

,  subcell  I  Tr  /  .  subcell \ 

$A  L 

V  7=1 


(58) 


where  V  is  the  parent  cell  volume,  Vi  is  the  subcell  volume,  and  ((|f“*c<2Z/).  is  the  cell 
average  flux  moment  for  the  i*  subcell.  The  subcell  exiting  flux  moment  is  assembled  as 
an  area  weighted  average  of  the  subface  values: 

..  #  subfaces 

i >r-7  e  4(*n  (59) 

A  ;  =  1 
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where  A  is  the  parent  exiting  face  area,  Aj  is  the  sub-face  area,  and  is  the 

zeroth  exiting  flux  moment  for  the  k  sub-face. 

3.6.2.  Linear  Characteristic  Method  (LC) 

The  LC  method  introduces  further  complexity  because  it  carries  first  moments  of 
each  distribution.  The  first  step  in  applying  the  LC  spatial  quadrature  is  to  determine  the 
coefficients  of  the  source  distribution.  A  convenient  coordinate  system  for  computing  the 
source  coefficients,  because  it  simplifies  the  form  of  the  moments  integrals,  is  shown  in 
Figure  5  (note  that  although  the  orientation  is  similar  to  the  local  case  1  coordinate 
system,  the  source  system  is  defined  for  the  parent  tetrahedron). 

The  source  moments  are  known  either  from  initial  conditions  or  from  having  been 
computed  from  the  previous  iteration  scalar  flux  values.  The  moments  must  be 
transformed  from  the  centroid  system  to  the  source  system  using 


(  c  ) 

■V 

K.j 

~XcO 

Sr' 

=  R 

- 

-yco 

SA 

kJ 

ki 

~  ZcOj 

where  R  is  the  rotation  matrix  constructed  in  the  same  way  as  equation  (43),  xc0,  yc0,  and 
zc0  are  the  coordinates  of  the  centroid  in  the  global  system,  and  x1?  y},  and  zx  are  the 
coordinates  of  the  source  system  origin  in  the  global  system  (since  all  nodes  are  defined  in 
the  global  system,  the  node  defining  the  source  system  origin  must  first  be  translated  to  the 
centroid  system).  In  the  source  system,  the  source  moments  are  defined  as 
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S(uJ ,v/ ,w') 


(61) 


source 


where  the  source  volume  integral  operator,  ( •  )source ,  has  the  same  form  as  the  subcell 
volume  integral  operator  given  in  equation  (46),  and  S(u' ,v',w  ')=a'+b'u '+c !v '+d'w ' . 
Evaluating  the  integrals  gives 


SA=a'  +  -jK  +«2  +W3)  +  j(V2  +v3>  +  ^~W3 


,/ 


S  /  =  —w  +  —  (u.  +  w,+2w3)w  +  — (v,+2v3)w3  +  —w2 
w  4  3  20  1  2  3  3  20  2  3  3  10  3 


d  ^  b  ^ 

Sv/  =  -j-(V2+Vl)  +  ^(UlV2+2u2V2+U3V2+UlV3+U2V3+2U3Vl) 
+  J^(V2  +V2V3+V3)  +  ^(V2+2V3>3 


(62) 


c  a'  ,  Z>7..  2  2  2s 

-  —  («J  +M2  +M3)  +  —  («j  +  UlU2  +W2  +Z/1W3  +  M2M3  +  U3  ) 


+  20  V2  +1U2V2  +  U3V2  +M1  V3  +  W2V3  +2«3V3) 


+  — (u,  +u ,  +2m,)w, 
20  1 
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Equations  (62)  define  a  system  of  four  equations  in  the  four  unknowns  a',  b',  c',  and  d' 
which  are  analytically  inverted  to  give  a  set  of  equations  that  is  used  to  directly  compute 
the  coefficients  given  the  source  moments. 

Once  computed  in  the  source  system,  the  coefficients  are  transformed  to  the 
centroid  system  using 


V 

Cc 

II 

cf 

A 

ff'i 

(63) 


and. 


a  =  a' -b't.-c't.-d't^ , 

C  1  2  3  7 


(64) 


where  R  is  the  rotation  matrix  used  to  transform  the  moments  to  the  source  system  and 


(  \ 
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h 

X1  ~XcO 

1 2 

a 

a 

» -yco 

As  j 

■  Z1  ~ZcO) 

(65) 


The  coefficients  in  the  centroid  system  can  then  be  transformed  into  the  local  sub- 
tetrahedra  systems  as  needed.  This  is  more  computationally  efficient  than  transforming  the 
source  moments  into  the  subcell  systems  and  directly  computing  the  coefficients  of  the 
distribution  for  each  subcell. 
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We  use  a  similar  approach  to  determine  the  coefficients  of  the  entering  flux 
distribution.  The  moments  are  known  in  the  global  face  system;  either  from  boundary 
conditions  or  from  upstream  cells.  In  the  global  face  system,  the  moments  of  the  entering 
flux  distribution  are  defined  as 


/  .  \ 

Va 

/ 

1 

1  ln 

♦x 

-\ 

Xp 

i  w 

Sr, 

r(VF> 


(66) 


where  t|f”  (xF,yF)  =  a  +  PxxF  +  (3 ~yyF,  and  the  entering  face  integral  operator,  ( •  )^ace ,  is 
defined  as 


(67) 


y  b  f  1  ~  Y ) 

where  x,  (yF)  =  — — yF  ,  and,  xR(yF)  =  b  -  — - —yF  •  Evaluating  the  integrals  gives 

h  h 


^“  =  i[3a+Z>px(1+Y)+/2py] 

=  j^[4a  +  *PX0  +2Y)  +  2*Py] 


4a(l+y)  +  2&PX(1  +  y  +y2)  +2/^(1  +  2y)  j  • 


(68) 
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Equations  (68)  comprise  a  linear  system  of  three  equations  in  the  three  unknowns  a,  px, 
and  Pv.  Solving  this  system  analytically  gives  a  system  of  equations  for  directly  computing 
the  coefficients  given  the  entering  flux  moments. 

Once  the  coefficients  of  the  source  and  entering  flux  coefficients  are  known,  they 
can  be  transformed  into  the  local  uvw  system  to  use  in  computing  the  subcell  moments. 
For  the  source  coefficients: 


V 

V 

c 

=  R 

C c 

\  d ) 

A 

(69) 


and  a  =  ac+btx+ct2  +dt3  -  bcxc0-ccyc0~dczc0.  The  entering  flux  coefficients  are 
transformed  from  the  global  face  system  to  the  local  uvw  system  using 


p, 

A, 


=  R 


face 


v 


(70) 


and  0!w  =  a  -  Pvtj  -  Pw/2 ,  where  in  this  case,  (^  ,/2)Tis  the  translation  vector  from  the 
global  face  to  the  local  uvw  origin  (the  translation  can  be  done  in  two  dimensions  because 
the  origins  of  both  systems  are  in  the  same  face). 

As  in  the  SC  method,  the  subcell  flux  moments  for  the  LC  method  are  the  sum  of 
the  source  and  entering  flux  contributions: 
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The  source  contributions  to  the  cell  flux  moments  are  defined  as 
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(72) 


where  S(s,  v,w)  =  a+bs+cv+dw  and  the  subcell-volume  integral  operator  was  defined 
in  equation  (46).  Evaluating  the  integrals  gives 
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The  entering  flux  contributions  to  the  cell  flux  moments  are  defined  as 


where  tyin[uin(v,w),v,w]  =  avw  +  pvv  +  py);w ,  yielding 
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t tin  ~  3  aw  M2  ( e )  +  [  Pv ( v2  +  v3 )  +  Pww3  j A/3  ( e ) 

=  «vw  W3M3  (  6  )  +  y  [  Pv(  V2  +  2  V3  )  +  2  PWW3  ] M4(e  ) 

tv"  =  «vw  (  V2  +  V3  )  M3(€  )  +  ^  [  Pv ( 2  V  +  2  V2  V3  +2  V32 )  +  PWW3  (  V2  +  2v3>]  (75) 

t  uj*11  =  3  Kvw  M1M2  ( 6  )  +  [  01  (  “3  M1  +  M2  +  «3  )  +  Pv"l  ( V2  +  V3  )  +  Pw^S  “l  ]M3  0 6  ) 

+  ^  [  Pv(2w2V2_4MlV2  +W3V2  ~4miV3  +W2V3  +  2m3V3) 

+  PwW3("4wi+W2+2W3)3M4(e)  • 
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where  the  face  integral  operator,  (m)face,  was  defined  in  equation  (54).  Evaluating  the 
integrals  gives  the  equations  for  the  source  contributions: 


The  entering  flux  contributions  to  the  exiting  flux  moments  are  defined  as 


which,  when  evaluated  give: 

CL  =  2  aM\  (  6  )  +  [Pv(v2  +  v3  )  +  }M2  (e ) 

CL  =  ^M2{e)  + y[Pv(v2  +2v3)  +  2§ww^M2{e) 

CL  =  «(v2  +  V3)M2(e)+l[Pv(2v22+2v2v3  +  2v32)  +  Pww3(v2+2v3)]M3(e)  . 
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The  subcell  moments  must  be  transformed  back  to  the  parent  cell  coordinate 
systems  before  they  can  be  assembled  into  the  parent  cell  results.  The  cell  flux  first 
moments  are  rotated  back  to  the  centroid  coordinate  system  using 
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where  R;  is  the  rotation  matrix  for  the  i*  sub-tetrahedron.  The  exiting  face  moments  are 
transformed  using 
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where  ( tx , t2)J  is  the  translation  vector  from  the  local  uvw  origin  to  the  face  system 
origin. 

In  the  centroid  system,  the  subcell  flux  moments  have  to  be  assembled  into  the 
parent  cell  results.  As  they  were  for  SC,  the  parent  cell  moments  are  just  the  volume 
weighted  average  of  the  subcell  moments: 
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where  V  is  the  parent  cell  volume,  Vj  is  the  subcell  volume,  and  ij;1  is  the  flux  moment  for 
the  i*  subcell.  Likewise,  the  exiting  flux  moments  for  the  parent  cell  are  just  an  area 
weighted  average  of  the  subcell  results: 
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(84) 


where  A  is  the  parent  exiting  face  area  and  \  is  the  subcell  face  area. 

3.6.3.  Exponential  Characteristic  Method  (EC) 

The  implementation  of  the  EC  method  is  similar  to  that  of  the  LC  method  in  that 
both  carry  first  moments  of  the  distribution  so  the  application  of  the  quadrature  includes 
the  same  sequence  of  steps.  There  are,  however,  several  significant  differences.  In  the 
equations  defining  the  moment  equations  for  the  LC  method,  the  only  exponential  term 
came  from  the  attenuation  of  the  entering  flux  and  source  distributions,  and  the 
integrations  produced  exponential  moment  functions  of  a  single  dimension.  The  EC 
distributions  are  defined  in  terms  of  exponentials,  so  the  equations  defining  the  moments 
generate  multi-dimensional  exponential  moment  functions.  Further,  the  systems  of 
equations  for  the  source  and  entering  flux  moments  are  nonlinear,  and  the  distribution 
coefficients  must  be  found  using  Newton's  method. 
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As  in  the  LC  method,  the  first  step  in  applying  the  EC  quadrature  is  to  transform 
the  source  moments  to  the  source  coordinate  system.  This  is  accomplished  using  equation 
(60).  In  the  source  coordinate  system,  the  moments  of  the  source  distribution  are  defined 


in  equation  (61)  where  S(u! ,v/ ,w')  =  a'  e 


b'u'+c'v'+d'w' ,  whjch  when  evaluated  gives 


SA  =  6a' eY  M0(X,Y,Z) 

Si  =  SA  w3  p(X,Y,Z) 

Si  =  SA  [v3  p (X,Y,Z)  +  v2  p (Y-X,Z-X,-X)] 

Si  =  SA  [»3  p(X,Y,Z)  +  u2  p(Y-X,Z-X,-X)  +  ux  p (Y-Z,X-Z,-Z)\ , 


(85) 


where 


Y  =  b'u3  +  c'v3  +  d'w3, 

X  =  Y  -  b'u2  -  c'v2,  (86) 

Z  -  Y  -  b'ux . 


As  in  the  LC  method,  the  system  defined  by  equations  (85)  can  be  inverted  to  find  the 
distribution  coefficients.  However,  unlike  LC,  this  system  is  non-linear  and  must  be 
solved  using  some  root  solving  method  such  as  Newton's  Method.  Appendix  B  contains  a 
complete  discussion  of  the  root-solving  procedure.  Once  computed  in  the  source  system, 
the  coefficients  are  rotated  back  to  the  centroid  system  using  equation  (63)  and 
a  =  a ' e  h  ?1  c  ‘2  d  ‘3  where  they  are  “stored”  until  needed  for  the  subcell  computations. 
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The  moments  of  the  entering  flux  distribution  for  the  EC  method  are  defined  in 


equation  (66)  where  i|r'”(xF, yF)  =  ae 


I \xr  +  p yVf 


.  Evaluating  the  integrals  gives 


C  =  2  ae^M0(Xin,YJ 
tf  =  rAh9{XinJin) 


C  =  € \yh  P(xin, YJ  +  b  p{xm - 4, -4) 


where  Xm=h$y  +  yb$x  and  Yin  =Xm~b§x. 


Equations  (87)  define  a  non-linear  system  in  the  constants  and  Ym  that  are 
solved  using  Newton's  method.  Appendix  B  contains  a  complete  discussion  of  the  two- 


dimensional  root  solving  procedure. 

Once  computed,  the  source  and  entering  flux  coefficients  are  transformed  to  the 
local  uvw  system.  For  the  source  moments,  the  first  moments  are  rotated  exactly  as  the 
LC  moments,  using  equation  (63).  The  translation  is  accounted  for  in  the  constant  term, 
which  in  this  case  is  a  =  a  e  bt'  *ct2+d‘3  bcX'°  c'y'°  d‘Zc0 _  The  entering  flux  first  moments  are 


rotated  using  equation  (70)  and  the  translation  is  accounted  for  in  a.  =  ae 


Pv'l  -Pwf2 


The  source  and  entering  flux  contributions  to  the  subcell  flux  moments  are  given 
by  equations  (72)  and  (74)  where  the  exponential  forms  of  the  distributions  are  used: 


■S'(s,v,w)  =  ae 


bs+cv+dw 


\\r,n[u.(v,w),v,w]  =  ctvwe 


PvV+PwW 


Evaluating  the  integrals  gives  the  equations  for  computing  the  source  contributions. 
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and  the  equations  for  computing  the  entering  flux  contributions. 
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The  equations  for  computing  the  source  and  entering  flux  contributions  to  the 
exiting  flux  moments  for  the  EC  method  are  generated  by  substituting  the  exponential 
forms  of  the  distributions,  equations  (32),  into  equations  (77)  and  (79).  This  gives  the 
source  contribution  equations. 


CL  -  2«1aey-M()(X,,ZJ.rj) 

CL  - 


CL  =  CLb  P<CCr,)  +  vjP(z,-.r,,r, 


(93) 


and  the  entering  flux  contribution  equations. 


CL  -  2aer-Mc(Xlm!,Z[J 

C»  =  CLb  (94) 

CL  *  <C[v,P v2P(Z,*-Cc.-C,)]  ■ 
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The  remaining  steps,  transforming  the  subcell  flux  moments  to  the  centroid  system, 
transforming  the  subcell  exiting  flux  moments  to  the  face  system,  assembling  the  parent 
cell  results,  and  the  equations  used  to  accomplish  them,  are  identically  the  same  for  both 
the  LC  and  EC  method. 

3.6.4.  Constant  Face  Variants 

Two  additional  methods  were  derived  and  tested:  Constant  Linear  Characteristic 
(CLC)  and  Constant  Exponential  Characteristic  (CEC).  These  methods  retain  the 
assumed  form  of  the  source  distribution  used  in  the  LC  and  EC  methods,  but  assume  a 
constant  distribution  on  the  face. 

3.6.4. 1.  Constant  Linear  Characteristic 

The  general  flow  of  the  CLC  method  is  the  same  as  that  for  the  methods  already 
presented.  First,  the  coefficients  of  the  source  moments  have  to  be  determined.  Both  the 
LC  and  CLC  methods  use  the  same  assumed  form  of  the  source  distribution,  so  the 
procedure  for  computing  the  CLC  source  distribution  coefficients  is  identical  to  the  LC 
method.  The  entering  flux  distribution  is  assumed  constant,  hence  it  is  handled  the  same 
way  as  the  SC  method. 

In  the  local  coordinate  system,  the  source  contributions  to  the  subcell  flux 
moments  are  computed  as  in  the  LC  method,  using  equations  (72).  The  equations  for 
computing  the  entering  flux  contributions  to  the  subcell  flux  are  generated  by  substituting 
the  constant  form  of  the  entering  flux  distribution  into  equation  (74)  giving 
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,  subcell  o  i  in  \  /  \ 

$A,in  =3^M2(e) 

.  subcell  .in  ns/  \ 

=^w3M3(e) 

(95) 

,  subcell  tin  /  \  n  s  /  \ 

tv,/»  =t.4(V2+V3)M3(e) 


subcell 


=3i|^i/1M2(€)  +^M3(e)(-3i/1+«2+M3)  . 


The  exiting  face  distribution  is  constant,  only  the  zeroth  moment  is  needed.  The 
source  contribution  to  the  exiting  flux  moment  is  computed  using  the  LC  source 
contribution  to  the  exiting  face  average  moment  in  equation  (78), 

u 

+  y[Z>(«i+M2+M3) +c(v2+v3)+Jw3]M3(e)  .  (96) 

The  entering  flux  contribution  to  the  exiting  flux  moment  is  the  same  as  for  the  SC 
method: 


(97) 


3.6.4.2.  Constant  Exponential  Characteristic 

The  source  moments  and  coefficient  calculations  are  the  same  for  the  EC  and  CEC 
methods.  The  entering  flux  moment  for  the  EC  method  is  handled  the  same  way  as  the  SC 
method.  In  the  local  coordinate  system,  the  source  contributions  to  the  cell  flux  moments 
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are  again  computed  using  equations  (89).  Because  the  entering  flux  is  assumed  constant, 
the  entering  flux  contributions  to  the  cell  flux  moments  are  the  same  as  for  CLC;  equations 
(95).  The  source  contribution  to  the  exiting  flux  moment  is  the  same  as  the  EC  source 
contribution  to  the  exiting  flux  average  moment  from  equation  (93), 


=  2^aer-M<s(Xs,Zs,V!),  (98) 

and  the  entering  flux  contribution  to  the  exiting  flux  moment  is  identical  to  that  used  in  the 
SC  and  CLC  methods. 

3.7.  Balance  Equations 

The  balance  equations  are  obtained  by  taking  moments  of  the  Boltzmann  transport 
equation  over  the  case  1  cell, 

(V 

u 
v 

where  the  volume  integration  is  carried  out  over  the  case  1  volume.  Using  equation  (98), 
the  zeroth-moment  equation  is 

3(^7  ~  O  +  eiT4  =  u1SA  ,  (100) 


3u 


\Jr(M,v,w)  +a(w,v,w)iji(w,v,vt')  =S(u,v,w ) 
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(99) 


cell 


the  w-moment  equation  is 


70 


3(C  -  C)  +  ef  =  U^' 


(101) 


and  the  v-moment  equation  is, 

3(C'  -  O  +  er  -  uj*.  (102) 


To  derive  the  u-Moment  equation  refer  to  Figure  1  and  recall  that  on  the  entering  face. 


z  x  in  in 

Ui„(V’W)=C  1  V  +  C2  W  > 


and  on  the  exiting  face. 


/  x  OUt  OUt 

Uout(V’W)=C  1  V  +  C2 


where 


w„ 


c\  ~ 


■>  C2  ~ 


V2w3 


°ut=U 2_Jh  md  out  = 


v2(w3-Wi)+v3(w1-m2) 


V2W3 


(103) 


(104) 


(105) 


The  u-Moment  equation  is  then 


^  /  out  ,  out  out  ,  ou\  ~  /  in  .  in  in  ,  /«x  ~  ,  out  .a  .  u  0  u  /i 

3(cl  i|fv  +c2  t|rw  )-3(Cj  i|fv  +c2  i|;w)  +  3m]^  -  +ei|r"  =  u1Su  .  (106) 


Using  equations  (103)  and  (104), 
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(107) 


.  in  in  ,  in  in  .  in 

=  C\  +  ^2  ^ 


out  ,  out 


out  .  out 


,  out  out  .  out  out  .  out  ,  out 

=  Cl  tv  +  tw  +  wit.4 


and  the  u-Moment  equation  can  now  be  written  in  a  form  consistent  with  the  other 
conservation  equations  as 


/  ,  out 

3(t„ 


O  -  w,  +  et" 


u lS1 


(108) 


The  balance,  or  conservation  equations  serve  several  purposes.  They  serve  as  a 
check  of  the  analytic  results  of  the  derivations.  The  analytic  forms  of  the  moment 
equations  for  each  of  the  quadratures  were  checked  by  substituting  them  into  the  balance 
equations.  During  execution,  after  the  cell  and  exiting  flux  moments  are  computed  for  a 
sub-tetrahedron,  the  results  are  checked  by  substituting  the  numerical  results  into  the 
balance  equations.  This  ensures  that  no  loss  of  digits  due  to  errors  or  numerical  ill- 
conditioning  corrupts  the  solution.  Alternatively,  the  equations  can  be  used  as  a  more 
efficient  means  of  computing  the  cell  flux  moments. 

Computing  the  cell  flux  moments  for  the  EC  method  using  equations  (89)  and  (91) 
is  very  expensive.  Using  the  conservation  equations  in  the  form 
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t.4 


«iSA  _  3JT 

e  e 

“A  +  3i|rw  _  3t|tw 
e  e 

^A  +  3i|fv  _  3i|fv 
e  e 

"A  +  fiA  +  3jC 

e 


~  ,  OK? 

3^„ 


(109) 


is  a  much  more  efficient  way  to  generate  the  cell  flux  moments,  since  all  of  the  quantities 
on  the  right  hand  side  of  the  equations  are  either  input  to  the  quadrature  (SA,SU,SV,SW), 
computed  as  part  of  the  quadrature  C'AA  or  would  have  been  computed 

anyway  to  check  conservation  on  the  subcell  ( i|r™ ,  ty™ ,  )  (since  only  the  coefficients  of 

the  entering  flux  distribution  are  transformed  to  the  subcell,  the  moments  over  the  subcell 
inflow  face  must  be  constructed  from  the  fitted  distribution).  Equations  (109)  become 
numerically  ill-conditioned  as  the  cell  optical  thickness  approaches  zero,  and  they  must  be 
used  with  care.  Further  details  on  the  use  of  the  conservation  equations  to  compute  the 
cell  moments  can  be  found  in  Appendix  D. 
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IV.  Arbitrary  Tetrahedra  Code  Algorithm 

The  characteristic  methods  were  implemented  in  the  TetSN  code  developed  as  part 
of  this  research.  The  TetSN  algorithm  is  a  direct  extension  of  Miller's  algorithm  (Miller, 
1993)  to  three-dimensions. 

4.1.  Initialization 

When  executed,  the  code  first  reads  runtime  data  from  a  file  containing  the  names 
of  the  input  and  output  files,  the  order  of  the  angular  quadrature,  the  spatial  quadrature  to 
use,  moments  function  computation  tolerance,  and  the  problem  convergence  tolerance. 

The  input  file  contains  all  information  pertaining  to  mesh  and  materials.  It  consists  of  a  list 
of  apex  triplets  (the  xyz  coordinate  of  each  apex  in  the  global  coordinate  system),  a  list  of 
node  numbers  associated  with  each  tetrahedron,  the  material  number  of  each  tetrahedron, 
the  material  cross  section  and  source  information  for  each  material  and  energy  group,  the 
group-to-group  scattering  cross  sections,  and  the  boundary  conditions.  The  code  then 
generates  lists  of  neighbor  tetrahedra,  boundary  tetrahedra,  boundary  faces,  face  areas,  cell 
volumes,  and  centroid  locations.  Given  the  angular  quadrature  order,  the  code  fills  the 
direction  cosine  and  weight  arrays  and  sets  the  number  of  angles.  The  highest  energy 
group  scalar  fluxes  are  then  initialized  either  to  an  a  priori  estimate  or  a  default  value  of 
zero. 

4.2.  Iteration 

Assuming  no  up  scatter  (particles  may  not  gain  energy  in  collisions),  the  group 
source  moments  are  computed  using  equation  (28).  Each  energy  group  uses  all  of  the 
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higher-energy-group-scalar  fluxes  to  compute  its  scattering  source.  Therefore,  the 
program  execution  consists  of  the  outer  iteration,  the  iteration  on  the  between  the  group 
scattering  source,  and  the  inner  iteration,  the  iteration  on  the  within  group  scattering 
source.  The  outer  iteration  proceeds  through  each  of  the  energy  groups,  beginning  with 
the  highest  energy  group.  The  scalar  flux  for  that  group  is  converged  (the  inner  iteration) 
and  then  used  to  update  the  down-scatter  source  for  each  of  the  lower-energy  groups. 

The  iteration  then  proceeds  down  through  the  rest  of  the  groups.  Note  that  for  down 
scatter  only,  the  outer  iteration  simply  consists  of  a  single  pass  through  each  of  the  energy 
groups.  If  up  scatter  were  allowed,  the  lower  group  fluxes  would  contribute  to  the 
scattering  source  for  higher  energy  groups,  and  the  outer  iteration  would  have  to  be 
repeated  until  the  group  fluxes  all  converged. 

For  the  first  order  methods,  the  coefficients  of  the  source  distribution  must  be 
computed  for  each  tetrahedron.  Because  the  scatter  is  assumed  to  be  isotropic,  the  source 
moments  are  independent  of  angle,  and  the  source  moment  computation  need  only  be 
performed  once  per  iteration  per  tetrahedron.  The  scalar  flux  moments  and  hence  the 
source  moments  are  computed  with  respect  to  the  centroid  in  the  centroid  coordinate 
system,  and  transformed  to  the  local  uvw  coordinate  system  where  the  coefficients  are 
determined.  Once  computed  in  the  source  coordinate  system,  either  by  a  linear  system  of 
equations  for  LC  and  CLC  or  by  Newton's  method  for  EC  and  CEC,  the  source 
distribution  coefficients  are  transformed  back  to  the  centroid  coordinate  system  where 
they  are  “stored”  until  needed. 
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At  the  start  of  each  inner  iteration,  the  code  saves  the  previous  iteration  scalar 
fluxes  to  use  to  check  for  convergence  at  the  end  of  the  current  iteration.  Before  each 
spatial  walk,  the  code  resets  all  cell  angular  fluxes  to  zero  and  sets  face  angular  fluxes  in 
accordance  with  the  boundary  conditions.  It  then  computes  the  dot  product  between  the 
streaming  direction  and  the  outward  normal  of  every  face  of  every  tetrahedron  to 
determine  whether  the  face  is  input  or  output.  A  push-down  stack  is  used  to  store  the 
indices  of  cells  for  which  all  necessary  input  data  are  available.  Every  tetrahedron  in  the 
mesh  is  compared  to  the  input  data  list  to  determine  which  tetrahedra  have  the  necessary 
input  information.  Tetrahedra  having  input  information  for  each  input  face  are  appended 
to  the  stack,  and  the  stack  pointer  is  incremented. 

After  all  the  tetrahedra  have  been  checked  for  input  data,  the  last  tetrahedron  is 
pulled  from  the  stack,  and  the  stack  pointer  decremented.  The  entering  flux  distribution 
coefficients  for  each  input  face  are  then  computed  using  the  equations  derived  in  Chapter 
3.  As  with  the  source  moments,  the  coefficients  of  the  entering  flux  distribution  must  be 
generated  from  knowledge  of  the  zeroth  and  first  moments.  All  face  flux  moments  are 
computed  in  terms  of  a  face  coordinate  system.  This  trivializes  the  passing  of  moments 
between  adjacent  tetrahedra.  Since  they  are  defined  in  the  same  coordinate  system,  the 
entering  face  moments  of  one  tetrahedron  are  identically  equal  to  the  exiting  flux  moments 
of  the  upstream  tetrahedron.  Once  computed,  the  entering  face  distribution  coefficients 
are  rotated  from  the  face  to  each  of  the  local  uvw  systems.  The  advantage  of  this  is  that 
the  coefficients  only  need  be  computed  (or  rootsolved)  once  for  each  input  face.  If  the 
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face  is  subdivided  among  sub-tetrahedra,  as  in  cases  2,  3,  and  4,  the  coefficients  are 
transformed  into  the  local  coordinate  system  for  each  surface. 

After  the  entering  face  coefficients  have  been  computed,  the  splitting  routine  is 
called  and  the  sub-tetrahedra  coordinate  systems  and  rotation  matrices  are  generated, 
fn  the  local  uvw  system,  the  cell  and  exiting  flux  moments  can  be  computed  using  the 
appropriate  spatial  quadrature  and  the  equations  derived  in  Chapter  3.  Once  computed, 
the  moments  of  the  cell  flux  distribution  for  each  sub-tetrahedron  are  transformed  into  the 
centroid  coordinate  system  where  they  are  combined  with  the  moments  from  other  sub- 
tetrahedra  to  form  the  parent  cell  flux  moments.  The  exiting  face  moments  for  each  sub¬ 
tetrahedron  are  transformed  into  the  global  face  system.  For  case  -2,  -3  and  case  4 
tetrahedra,  the  exiting  face  moments  for  several  tetrahedra  must  be  combined  into  the 
parent  face  results. 

After  the  parent  tetrahedron  results  are  assembled,  the  current  angle  contributions 
to  the  scalar  flux  and  currents  for  the  tetrahedron  are  added  to  running  totals.  The  exiting 
face  flux  moments  are  used  to  update  the  input  information  for  neighbor  tetrahedra.  The 
neighbor  tetrahedra  are  then  checked  to  see  if  they  now  have  the  necessary  input 
information  and,  if  so,  are  appended  to  the  stack  and  the  pointer  incremented.  This 
process  is  repeated  until  the  stack  pointer  decrements  to  zero,  indicating  that  the  stack  is 
empty,  and  all  spatial  cells  have  been  computed  for  the  current  angle.  A  new  angle  is  then 
chosen  and  the  process  repeated  until  all  angles  in  the  quadrature  have  been  used.  The 
current  iteration  scalar  flux  is  compared  to  the  previous  iteration  scalar  flux  for  each 
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tetrahedron  to  determine  if  the  solution  has  converged.  The  comparison  for  each 
tetrahedron  is  performed  by  computing  an  error  ratio,  eR,  defined  as 


«*  =  2 


1 4>r  -  $  i 
i4>n+i4>ii  ’ 


(110) 


where  the  superscript  k  indicates  the  iteration  number.  If  the  largest  error  ratio  is  smaller 
than  the  specified  convergence  tolerance,  the  inner  iteration  is  converged,  otherwise  the 
current  group  source  is  updated  using  the  just-computed  scalar  fluxes,  and  the  process 
repeated  until  the  group  scalar  flux  (within  group  scattering  source)  converges.  When  the 
inner  iteration  converges,  the  current  group  scalar  flux  is  used  to  update  the  down  scatter 
source  for  the  lower  energy  groups.  A  new  inner  iteration  begins  using  the  next  lowest 
energy  group  and  continues  until  all  of  the  energy  groups  have  been  used.  At  this  point, 
execution  is  complete  and  the  output  file  is  generated. 

C.  Pseudocode  Algorithm 

The  following  is  a  Pseudocode  outline  of  the  TetSN  algorithm: 

Read  input  data 
Initialize  angular  quadrature 
Compute  face  normals 

Initialize  scalar  flux  arrays  for  current  problem 
For  each  energy  group: 

Do 
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Save  previous  iteration  scalar  fluxes,  reset  scalar  flux  arrays 

Update  the  source  moments  for  current  group 

Generate  source  coefficients  for  each  tetrahedron  if  LC, 

EC,  CEC,  orCLC 

For  each  angle: 

Reset  angular  flux  moment  arrays  to  zero 

Compute  dot  products  of  face  normals  and  current 
streaming  direction 

Determine  which  tets  are  ready  to  compute  and 
append  to  stack 

While  stack  not  empty: 

Take  tetrahedron  from  stack,  decrement 
pointer 

Determine  the  case  of  the  current 
tetrahedron 

Split  along  streaming  direction  and  generate 
rotation/translation  matrices  for  each  sub¬ 
tetrahedron 

Solve  for  the  coefficients  of  the  entering  flux 
distribution  for  each  input  face  of  the  parent 
tetrahedron 

For  each  sub-tetrahedron: 

Rotate  the  coefficients  of  the  source 
distribution  to  the  local  uvw 
coordinate  system 

Rotate  the  coefficients  of  the  input 
flux  distribution  to  the  local  uvw 
system 
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Compute  the  cell  flux  and  exiting 
flux  distribution  moments 

Rotate\translate  the  cell  flux 
moments  to  the  centroid  coordinate 
system 

Rotate/translate  the  exiting  flux 
moments  to  the  face  system 

Next  sub-tetrahedron 

Assemble  parent  cell flux  moments  from 
sub-tetrahedron  moments 

Assemble  parent  exiting  face  moments  from 
sub-face  moments 

Update  neighbor  tetrahedra  face  flux  values 
and  append  to  stack  those  now  ready  to 
compute 

Accumulate  scalar  flux  moments  and 
currents  for  current  tetrahedron 

Next  tetrahedron  from  stack 

Next  Angle 

Check  for  within  group  convergence  by  comparing 
previous  iteration  scalar  flux  to  current  iteration  fluxes 

Loop  until  Converged 

Compute  current  group  contribution  to  lower  energy  group  down  scatter 
sources 

Next  Energy  group 
Generate  output 
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V.  Testing 

The  characteristic  quadratures  in  the  TetSN  code  were  evaluated  on  four  test 
problems.  The  first  two  test  problems  consist  of  simple  geometries  that  by  no  means  test 
the  complex  three-dimensional  solution  capabilities,  but  do  demonstrate  basic  operation 
and  convergence  of  the  methods  without  unnecessary  clutter.  The  geometry  for  test 
problems  3  and  4  is  fully  three-dimensional.  Although  all  of  the  methods  are  tested,  heavy 
emphasis  is  placed  on  the  first-moment  methods,  LC  and  EC.  The  constant  face  methods 
operate  correctly,  but  their  low  order  of  convergence  and  lack  of  accuracy  limits  their  use 
in  practice,  and  as  a  result,  they  are  of  academic  interest  only. 

Test  problems  1  and  3  are  benchmarked  against  MCNP,  (Monte  Carlo  Neutron 
Photon  Transport  Code  System)  (Briesmeister:  1991).  The  MCNP  code,  developed  at 
Los  Alamos  National  Laboratory,  is  a  Monte  Carlo  simulation  code  for  neutral  particle 
transport  in  general  three-dimensional  geometry.  It  is  widely  used  and  accepted  in  the 
transport  community.  The  MCNP  result  is  the  average  of  the  computed  quantity  over  all 
histories,  x,  defined  as 


x  = 


(111) 


where  N  is  the  number  of  histories  calculated  in  the  problem  and  Xj  is  the  value  of 
computed  quantity  for  the  i1*1  history.  MCNP  also  reports  the  estimated  statistical  relative 
error  at  the  la  level,  R,  defined  as  R  =  S~/x  where  S2  is  the  estimated  sample  variance. 
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Test  problem  1  demonstrates  the  basic  operation  of  the  code,  and  correctness  of 
the  individual  methods.  The  convergence  rate,  sensitivity  to  mesh  variations,  and 
execution  speed  of  each  method  are  evaluated  on  test  problem  1 .  The  test  problem  1 
results  were  benchmarked  against  MCNP  solutions.  MCNP  was  used  as  an  independent 
check  to  ensure  that  there  were  no  gross  errors  in  the  derivations  and  implementations. 
Test  problem  2  is  a  deep  shielding  problem  used  to  demonstrate  the  superior  performance 
of  the  EC  method  on  very  thick  cells.  The  test  problem  2  results  are  benchmarked  using  a 
semi-analytic  solution  that  treats  the  spatial  variable  exactly  but  uses  a  discrete  ordinates 
angular  approximation.  Test  problem  2  has  slab  symmetry  in  order  to  make  the  semi- 
analytic  solution  feasible.  Test  problem  3  utilizes  true  three-dimensional  geometry  to 
demonstrate  the  effectiveness  of  the  methods  in  three  dimensions,  and  to  illustrate  the 
need  to  conserve  volumes  when  meshing  regions,  a  limitation  of  the  current  mesh 
generation  method.  Test  problem  3  results  are  benchmarked  against  MCNP  solutions. 
Test  problem  4  uses  a  three-dimensional  problem  to  illustrate  the  effectiveness  of  the  EC 
method  on  thick  cells  and  the  ability  to  use  the  code  to  generate  flux  maps  on  a  problem 
that  would  be  difficult  to  solve  using  Monte  Carlo.  Because  the  problem  contains  curved 
surfaces,  it  also  demonstrates  the  power  of  tetrahedral  meshes  over  traditional  boxoid 
meshes. 
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5.1.  Test  Problem  1  -  Unit  Source  Cube  with  Vacuum  Boundaries 


The  first  test  problem  is  a  simple  lxlxl  cm  cube  with  vacuum  boundaries  and  a 
uniformly-embedded  isotropically-emitting  source  of  strength  1.0  cm'J  sec'1.  The  cube 
material  has  a  scattering  cross  section  of  0.5  cm'1  and  a  total  cross  section  of  1.0  cm'1. 
The  meshes  for  test  problem  1  were  originally  generated  using  the  IDEAS  code,  but  the 
convergence  data  generated  using  them  was  contaminated  by  the  mesh  sensitivity 
variations  to  be  discussed  in  Section  5.1.2.  To  isolate  the  convergence  data  from  the 
sensitivity  variations,  a  series  of  structured  meshes  was  generated  manually.  To  construct 
each  mesh,  the  problem  domain  was  first  meshed  using  cubic  cells,  then  each  cube  was 
subdivided  into  six  tetrahedra.  This  resulted  in  a  series  of  structured,  self-similar  meshes, 
the  first  three  of  which  are  shown  in  Figure  6. 


Table  1.  Mesh  Sizes  Used  for  Test  Problem  1. 


Mesh  Number 

Number  of  Tetrahedra 

1 

6 

2 

162 

3 

1296 

4 

4374 

5 

10368 

6 

20250 

In  analogy  to  the  rectangular  case,  1  / («teK)1/3  is  an  approximation  to  Ax,  the 
average  linear  cell  dimension.  For  the  test  problem  1  mesh  sizes,  given  in  Table  1,  the 
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average  cell  thickness  in  mean  free  paths  ranges  from  0.55  to  0.036.  The  maximum 
optical  thickness  ranges  from  1.7  mfp  for  the  coarsest  mesh  to  about  0.1  on  the  finest 
mesh.  The  test  problem  1  solutions  were  computed  using  the  S16  level  symmetric 
quadrature  and  a  convergence  tolerance  of  10"6  . 


Figure  6.  First  Three  Test  Problem  1  Mesh  Refinements. 
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5.1.1.  Convergence 

There  are  two  levels  of  convergence  to  consider  in  evaluating  the  performance  of 
the  TetSN  code.  The  first  is  convergence  of  the  TetSN  solution  to  the  “correct”  solution, 
i.e.  the  result  that  would  be  obtained  in  the  laboratory.  Convergence  in  this  sense  is 
obviously  important  because  to  be  useful,  the  discrete  ordinates  solution  must  accurately 
approximate  reality.  Because  the  methods  use  discretized  spatial  and  angular  variables,  the 
solutions  generated  for  a  given  mesh  and  angular  quadrature  order  are  approximate.  In 
theory,  as  the  spatial  and  angular  domains  are  refined  the  discrete  ordinates  solution 
should  converge  upon  the  exact  solution. 

The  second  level  of  convergence  is  a  subset  of  the  first,  and  illustrates  the 
efficiency  of  the  spatial  quadrature  used  to  obtain  the  discrete  ordinates  solution.  In  this 
case  we  are  concerned  with  how,  for  a  given  angular  quadrature  order,  a  particular  spatial 
quadrature  approaches  the  approximate  numerical  solution.  Even  if  the  spatial  variable 
were  treated  exactly,  the  solution  would  only  be  an  approximation  due  to  the  angular 
discretization.  As  the  spatial  mesh  is  refined  while  the  angular  quadrature  order  is  held 
constant,  the  solutions  approach  a  semi-analytic  result  in  which  the  angular  variable  is 
discretized  and  the  spatial  variable  is  treated  exactly. 

As  the  mesh  is  refined,  the  average  linear  cell  length,  Ax,  decreases  by  a  factor 
n  =  Ax,  /A x2  where  Ax,  >  Ax2  are  the  average  cell  lengths  for  each  refinement.  In  the 
thin  cell  limit,  a  factor  n  reduction  in  the  linear  dimensions  of  a  spatial  cell  should  result  in 
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a  reduction  in  the  solution  error  given  by 


(112) 


where  €j  and  e2  are  the  solution  errors,  relative  to  the  benchmark  result,  for  the  two 
meshes,  and  p  is  the  order  of  convergence.  For  methods  with  a  large  order  of 
convergence,  on  fine  meshes,  further  refinement  rapidly  decreases  the  truncation  error  so 
that  greater  accuracy  is  easily  achieved. 

We  first  examine  the  convergence  of  the  spatial  quadratures  to  the  fully  spatially- 
converged  S16  discrete  ordinates  solution.  The  problem  average  flux  as  computed  by 
TetSN  using  each  of  the  quadratures  is  plotted  in  Figure  7  and  Figure  8.  The  exiting  face 
current  results  for  a  single  face  are  shown  in  Figure  9  and  Figure  10.  Because  of  the 
symmetry  of  the  problem,  all  faces  are  equivalent.  The  first  order  methods  appear  to  have 
converged  to  the  same  solution  on  the  finest  (20250  cell)  mesh.  The  constant  face 
methods  are  approaching  the  first  order  solution,  but  much  more  slowly.  This  is  better 
seen  by  examining  the  relative  error  at  each  mesh  refinement. 
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lu 


Figure  7.  Problem  Average  Flux  for  Test  Problem  1. 


(Enlargement  of  fine  mesh  results) 
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Figure  9.  Exiting  Face  Flux  for  Test  Problem  1. 


Figure  10.  Exiting  Face  Current  for  Test  Problem  1. 
(Enlargement  of  Fine  Mesh  Results) 
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In  Figure  11  and  Figure  12,  the  absolute  relative  error,  e^,  of  the  solution  for  each 
spatial  refinement  with  respect  to  the  benchmark  discrete  ordinates  solution  is  plotted 
against  the  cube  root  of  the  number  of  tetrahedra  in  the  mesh  for  each  of  the  methods. 

The  absolute  relative  error  between  a  given  value,  4>,  and  a  benchmark  result,  (f>bench ,  is 
computed  using 


_  |  §bench  -  (j) 

ar  Kbench 


<j> 


(113) 


Because  we  are  concerned  with  the  convergence  of  the  spatial  quadrature  to  the  S16 
discrete  ordinates  solution  as  the  mesh  is  refined,  we  use  the  most  converged  TetSN 
result,  the  20250  tetrahedron  mesh  LC  solution,  as  the  benchmark.  On  the  log-log  plot, 
the  order  of  convergence  for  a  method  is  simply  the  slope  of  its  corresponding  curve.  EC 
and  LC  both  exhibit  approximately  third  order  convergence.  This  is  in  agreement  with  the 
third-order  convergence  rates  observed  for  these  methods  in  two-dimensional  rectangular 
and  triangular  schemes  (Minor,  1993:  91-94;  and  Miller,  1993:74-76).  The  constant  face 
methods  however,  exhibit  slightly  less  (p~0.9)  than  the  first-order  convergence  rates 
observed  for  the  SC  method  (Minor,  1993:  91-94;  and  Mathews,  Miller  and  Brennan, 
1996).  As  discussed  in  Section  2.2  ,  using  constants  to  represent  the  entering  flux 
distribution  can  cause  numerical  diffusion.  The  constant  face  quadratures,  as  derived  for 
tetrahedral  cells,  appear  to  be  more  sensitive  to  the  treatment  of  the  entering  face  flux  (see 
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also  Section  5.1.2)  than  is  the  case  in  one  or  two  dimensions,  contributing  to  their  slow 
convergence  and  lack  of  accuracy. 

It  should  also  be  noted  that  the  MCNP  results  for  the  problem  average  flux  and 
exiting  current  are  inconsistent.  The  MCNP  result  gives  a  higher  problem  average  flux 
than  the  TetSN  solutions.  This  being  the  case,  we  expect  the  MCNP  exiting  current  to  be 
lower.  This  is  not  the  case.  Because  of  the  nature  of  the  estimators  used  in  the  Monte 
Carlo  calculation,  particle  balance  is  not  strictly  satisfied.  If  we  recompute  the  value  of 
the  exiting  current  based  on  the  MCNP  average  flux  result  and  conservation  of  particles, 
we  get  the  result  shown  in  the  dashed  line  of  Figure  9  and  Figure  10.  Indeed,  the  TetsN 
solutions  for  the  exiting  current  are  converging  toward  the  particle  conserving  MCNP 
results. 


Figure  11.  Relative  Error  in  Test  Problem  1  Problem  Average  Flux. 
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Figure  12.  Relative  Error  in  Test  Problem  1  Exiting  Face  Current. 


Error  norms  of  the  face  flux  distribution  were  computed  to  estimate  the  pointwise 
convergence  properties  of  the  quadratures.  The  error  norms  were  computed  by  mapping  a 
regular  mesh  of  points  onto  the  surface.  For  each  point,  the  flux  value  was  computed 
based  on  the  computed  flux  representation  for  the  tetrahedron  face  containing  the  point. 
This  allowed  us  to  compare  the  distributions  of  two  separate  meshes  having  dissimilar  face 
arrangements. 

Two  norms  were  computed:  and  L} .  The  norm  is  defined  as  the  maximum 

of  the  absolute  value  of  the  relative  errors  over  the  face, 


=  max 

n  =  ljf 


(114) 
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where  N  is  the  number  of  points  used,  ({)(x,  y)  is  the  flux  distribution  at  point  (x,  y)  for  the 
current  mesh  refinement,  and  ({^^(x,  y)  is  the  flux  at  point  (x,  y)  for  the  benchmark 
solution.  The  Lj  norm  is  defined  as  the  average  of  the  absolute  relative  errors  over  the 
face, 

_  i  ft  \<b(x„,y„)-<bbench(xn,yn)\ 

1  N  „=i  $bench(x„,yn) 

The  norm  computations  were  repeated  using  successively  finer  point  spacings  until  the  L2 
norm  converged  to  three  digits  and  the  L„  norm  converged  to  two  digits  (the  L„  norm  is 
more  sensitive  to  small  changes  in  the  point  spacing).  This  required  using  a  50x50  mesh 
of  points. 

Figure  13  shows  the  L„  norm  for  the  exiting  face  flux  distribution  compared  to  the 
benchmark  (20250  tetrahedron  LC  result)  at  each  spatial  refinement  for  each  of  the 
quadratures.  The  average  relative  error  of  the  exiting  flux  is  shown  in  Figure  14.  The 
first  order  method  norms  show  approximately  second  order  convergence.  Although  the 
constant  face  methods  do  a  reasonable  job  of  approximating  the  average  value  of  the 
exiting  flux,  we  should  not  expect  them  to  accurately  approximate  the  pointwise 
distribution.  The  constant  solution  on  a  face  agrees  with  the  pointwise  values  only  where 
the  pointwise  solution  equals  the  average.  Even  for  a  relatively  flat  distribution,  the 
constant  face  methods  require  a  fine  mesh  to  adequately  resolve  the  face  flux  spatial 
distribution. 
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Figure  13.  L-Infinity  Norm  (Maximum  Relative  Error)  of  Test  Problem  1 
Exiting  Face  Flux  Distribution. 


Figure  14.  L-l  Norm  (Average  Relative  Error)  of  Test  Problem  1 
Exiting  Face  Flux  Distribution. 
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The  Monte  Carlo  solution  method  uses  continuous  angular  and  spatial  variables 
and  contains  only  statistical  errors.  Because  of  the  relative  simplicity  of  the  problem 
geometry  and  cross  sections,  the  MCNP  solution  represents,  to  within  statistical 
uncertainty,  the  laboratory  solution.  The  MCNP  solution  is  also  plotted  in  Figure  7 
through  Figure  10.  The  uncertainty  in  the  MCNP  solution  for  the  average  flux  is  0.02%, 
which  is  too  small  to  represent  with  error  bars  in  Figure  7,  but  is  shown  in  Figure  8.  The 
MCNP  exiting  current  uncertainty  was  0.08%.  The  difference  in  the  fully  converged 
TetSN  and  MCNP  solutions  is  a  result  of  the  discretized  angular  treatment  used  in  the 
discrete  ordinates  calculation.  Overall,  we  expect  that  as  the  solution  approaches  the 
continuous  solution,  i.e.  as  the  discretized  spatial  and  angular  variable  mesh  spacing 
approaches  zero,  the  discrete  ordinates  solution  should  approach  the  actual  (in  this  case 
MCNP)  solution.  Figure  15  shows  the  relative  error  between  TetSN  EC  solutions  for 
various  mesh  refinements  and  angular  quadrature  orders  and  the  MCNP  solution.  Figure 
16  and  Figure  17  show  the  same  for  the  LC  and  SC  methods  respectively.  The  figures 
show  that  each  method  converges  to  the  MCNP  solution  as  either  the  angular  quadrature 
order  or  number  of  tetrahedra  in  the  mesh  is  increased,  as  it  should  for  a  properly  running 
code. 
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Figure  15.  Relative  Error  Between  TetSN  EC  and  MCNP  Solutions  for 
Test  Problem  1  Problem  Average  Flux. 


Figure  16.  Relative  Error  Between  TetSN  LC  and  MCNP  Solutions  for 
Test  Problem  1  Problem  Average  Flux. 
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Figure  17.  Relative  Error  Between  TetSN  SC  and  MCNP  Solutions 
for  Test  Problem  1  Problem  Average  Flux. 


5.1.2.  Mesh  Sensitivity 

A  robust  method  should  produce  consistent  results  on  a  mesh  containing  a  given 
number  of  cells,  regardless  of  the  particular  arrangement  of  the  individual  cells.  If  this  is 
not  the  case,  then  it  is  impossible  to  have  any  confidence  in  any  one  solution  obtained  on 
an  arbitrarily  generated  mesh.  To  study  the  effects  of  mesh  variations,  each  method  was 
used  to  generate  solutions  for  test  problem  1  on  a  series  of  ten  different  meshes,  each 
containing  162  cells.  The  modified  meshes  were  generated  by  randomly  varying  the  x,  y. 
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and  z  coordinates  of  each  node  in  the  original  mesh  (comer  nodes  were  fixed  and  the  face 
and  edge  nodes  were  constrained  to  variations  in  the  face  plane  and  edge  line  respectively 
so  as  to  retain  the  original  problem  boundaries).  The  node  coordinates  were  allowed  to 
vary  by  up  to  25%  or  50%  of  the  distance  between  the  two  closest  nodes  in  the  mesh,  i.e. 

^  =  ev.r^min(2v-1)+^g 

y™  =  evaA min  (2  V  "  1  )  +>W  (1  16) 

ZneW  =  €vardmm(2V  -l  )+Zong’ 

where  dmin  is  the  minimum  mesh  point  spacing,  evar  is  the  maximum  fractional  variation 
allowed  (0.25  or  0.5),  v  is  a  random  number  that  varies  between  0  and  1,  and  xorig,  yorig, 
and  zorig  are  the  original  node  coordinates.  This  produced  the  two  sets  of  meshes  shown 
in  Figure  18  and  Figure  19. 
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Figure  18.  Exterior  View  of  25%  Figure  19.  Exterior  View  of  50% 

Variation  Meshes.  (Apparant  distortion  Variation  Meshes, 
of  the  cube  is  artifact  of  file  transfer) 


Before  examining  the  influence  of  the  mesh  variations,  we  examine  the 
characteristics  of  the  meshes  themselves.  As  a  quantifiable  measure  of  tetrahedron  shape, 
we  define  the  aspect  ratio  of  a  given  tetrahedron  to  be  the  ratio  of  the  tetrahedron  volume 
to  that  of  an  equilateral  tetrahedron  in  which  the  edges  are  the  length  of  the  longest  edge 
of  the  tetrahedron  in  question.  For  this  metric,  an  ideal  (equilateral)  tetrahedron  has  a 
volume  ratio  of  one,  and  a  poorly  shaped  tetrahedron  has  a  volume  ratio  approaching 
zero.  The  cumulative  distribution  functions  (CDF)  of  volume  ratios  for  the  two  sets  of 
meshes  are  shown  in  Figure  20  and  Figure  21 .  The  Heaviside  distribution  in  the  figures  is 
the  structured  mesh  CDF,  for  which  all  tetrahedra  have  a  volume  ratio  equal  to  0.27. 


_ VolumeRatio _ 

Figure  20.  Cumulative  Distribution  Function  of  Tetrahedron  Volume  Ratios 
for  the  25%  Variation  Meshes. 
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Figure  2 1 .  Cumulative  Distribution  Function  of  Tetrahedron  Volume  Ratios 
for  the  50%  Variation  Meshes. 


The  effect  of  the  variations  was  to  produce  tetrahedra  approximately  equally 
distributed  about  the  structured  mesh  value.  The  50%  cases  contain  more  extreme 
variations,  including  some  very  poorly  shaped  tetrahedra.  For  the  ten  50%  meshes,  the 
minimum  volume  ratios  ranged  from  0.00017  to  0.019,  and  the  maximum  ranged  from 
0.53  to  0.9.  Figure  22  shows  the  structured  and  two  50%  variation  cases  (mesh  #4  and 
mesh  #5)  along  with  the  worst  and  best  shaped  tetrahedron  from  each.  The  tetrahedron 
shapes  ranged  from  essentially  plane  figures  to  almost  perfect  equilateral  tetrahedra. 

To  be  robust,  the  code  must  be  able  to  handle  wide  variations  in  the  cell  aspect 
ratio.  One  would  hope  that  a  mesh-generating  algorithm  would  at  least  avoid  such 
dramatically  bad  shapes.  This  is  indeed  the  case  for  the  IDEAS  generated  meshes.  A 
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Figure  22.  Exterior  Mesh  View  of  Structured  Mesh,  Mesh  #4  and  Mesh  #5  with 
Worst  and  Best  Shaped  Tetrahedron  for  Each. 

CDF  for  a  482  cell  mesh  generated  by  IDEAS  is  plotted  along  with  the  structured  mesh 
CDF  and  . a  representative  sample  from  each  of  the  25  and  50%  variations  in  Figure  22. 
The  tetrahedra  generated  by  IDEAS  have,  on  average,  better  aspect  ratios  than  the 
tetrahedra  generated  by  randomization.  Hence  for  most  problems,  IDEAS  should  be  able 
to  produce  acceptably  shaped  tetrahedra  and  there  is  no  need  to  be  concerned  with 
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method  performance  on  poorly  shaped  cells.  There  are,  however,  problems  for  which  the 
shapes  and  relative  sizes  of  the  individual  regions  in  a  problem  require  either  many  high 
aspect  ratio  cells,  or  a  few  low  aspect  ratio  cells,  to  model.  For  these  problems, 


Figure  23.  CDF  for  IDEAS  Generated  Mesh. 


algorithms  capable  of  producing  accurate  results  on  poorly  shaped  tetrahedra  offer 
increased  computational  efficiency. 

Figure  24  and  Figure  25  show  the  relative  error  in  the  test  problem  1  problem 
average  flux  for  each  of  the  50%  mesh  variations.  Variation  0  is  the  structured  mesh.  The 
same  is  shown  for  the  face  flux  error  norms  in  Figure  26  and  Figure  27.  The  25% 
variations  for  all  computed  quantities  were  smaller  than  the  50%  values,  so  we  examine 


102 


only  the  50%  variation  meshes.  While  the  constant  face  methods  show  some  variations, 
the  relative  errors  are  all  of  the  same  order  of  magnitude  or  less  than  the  structured  mesh 
solution.  The  first  moment  methods  were  less  sensitive  to  variations  in  the  mesh,  with  LC 
performing  slightly  better  than  EC.  The  first  moment  method  solutions  for  the  average 
flux  were  accurate  to  three  digits,  but  the  variations  produced  by  the  different  meshes 
caused  the  solutions  to  change  in  the  fourth  digit,  hence  the  effects  of  the  variations  for 
the  LC  and  EC  quadratures  were  entirely  insignificant.  Note  that  the  presence  of  the  very 
poorly  shaped  tets  in  meshes  4  and  5  caused  no  significant  deviations  in  the  relative  errors 
for  those  meshes. 

As  explained  in  section  2.2,  using  constants  to  represent  the  face  flux  distributions 
can  lead  to  truncation  errors,  namely  numerical  diffusion.  The  way  the  errors  generated  as 
a  result  of  numerical  diffusion  propagate  through  the  mesh  depends  on  the  specific 
orientation  of  the  faces,  and  hence  the  cells,  in  the  mesh.  This  means  that  the  constant 
face  method  errors  (more  so  than  the  first  moment  method  errors),  relative  to  the  most 
converged  discrete  ordinates  solution,  depend  upon  both  the  number  of  cells  and  the 
specific  orientation  of  the  cells.  This  heightened  sensitivity  to  variations  in  the  mesh 
composition  (Figure  24)  required  the  use  of  structured  meshes  to  evaluate  quadrature 
convergence  in  Section  5.1.1. 
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Figure  24.  Relative  Error  in  Problem  Average  Flux  for  50%  Node  Variation. 
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Figure  25.  Relative  Error  in  Problem  Average  Flux  for  50%  Node  Variation. 
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e la  t  i  v 


Figure  26.  L,  Norm  of  Exiting  Face  Flux  for  50%  Variation. 


Figure  27.  LI  Norm  of  Exiting  Face  Flux  for  50%  Variation. 
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5.1.3.  Execution  Speed 

The  execution  times  per  phase  space  cell  (i.e.  per  tetrahedron,  per  angle,  per 
energy  group,  per  iteration)  for  test  problem  1  on  an  IBM  RSC-6000,  model  590 
computer  are  given  in  Table  2. 


Table  2.  Execution  Time  per  Phase  Space  Cell  (seconds). 


Total 

Quadrature 

/Cnsrv. 

Split/ 

Assemble 

Source 

Moments 

Input 

Moments 

Over¬ 

head3 

SC 

1.2X10-4 

6.4x1  O'5 

3.5xl0'5 

- 

- 

2.  lxl  O'5 

CLC 

1.4X10-4 

7.3x1  O'5 

4.2x1 0'5 

1.7x1  O'7 

- 

2.5xl0'5 

CECb 

l.lxlO'3 

9.7x1  O'4 

5.0xl0"5 

2.9x1  O'5 

- 

2.9x1  O'5 

LC 

2.0x10^ 

7.7x1 0"5 

7.3x1  O'5 

1.7x1  O'7 

4.9x1  O’6 

3.5xl0"5 

ECC 

1.6xl0‘3 

1.0x1  O'3 

8.9xl0'5 

2.9x1  O'5 

3.6X10"4 

4.0x1 0'5 

a.  overhead  includes  reading  input;  generating  neighbor,  face,  and  centroid 
lists;  accumulating  scalar  fluxes;  computing  dot  products;  resetting  arrays 
between  iterations;  updating  neighbor  triangles;  and  writing  the  output 
file. 

b.  the  CEC  method  was  never  modified  to  use  the  balance  equations  to 
compute  the  cell  source  moments. 

c.  EC  (and  to  a  smaller  extent  CEC)  execution  times  may  vary  by  about 
±30%  from  the  listed  times  depending  on  the  details  of  the  specific 
problem  being  solved  (see  text) 

The  overhead  tasks  can  be  divided  into  four  categories:  those  required  to  be 
performed  once  per  program  execution,  those  required  once  per  iteration,  those  required 
once  per  angle  per  iteration,  and  those  required  once  per  cell  per  angle  per  iteration.  The 
overhead  tasks  required  once  per  cell  per  angle  per  iteration  account  for  about  86%  of  the 
overhead  time  in  the  problem  used  to  generate  the  table  (162  tets,  9  iterations,  80  angles). 
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Further,  the  majority  of  the  remaining  tasks  scale  with  the  number  of  spatial  cells.  Hence, 
the  time  required  to  execute  the  overhead  tasks  scales  with  the  number  of  cells ,  and  the 
times  in  the  table  can  be  considered  to  be  independent  of  problem  size. 

The  EC  method  is  the  most  costly,  requiring  about  8  times  the  effort  of  LC  and  14 
times  the  effort  of  SC  per  phase  space  cell  (per  tetrahedron  per  angle  per  iteration).  The 
main  difference  between  the  exponential  methods  and  the  other  methods  is  that  the 
exponentials  in  the  assumed  form  of  the  flux  and  source  distributions  produce  multi¬ 
dimensional  moment  functions  in  the  resulting  equations.  The  multi-dimensional  functions 
are  generally  more  expensive  to  compute.  The  SC,  CLC,  and  LC  quadrature  moment 
equations  all  contain  moments  functions  of  a  single  argument,  in  fact,  always  the  same 
argument,  e,  making  them  much  more  efficient.  The  EC  method  spends  63%  of  its 
execution  time  computing  moments  functions,  while  the  LC,  CLC,  and  SC  methods  spend 
between  20  and  30%.  Conversely,  the  non-exponential  quadratures  spend  53%  of  their 
execution  time  splitting  and  assembling  tetrahedra  while  the  EC  method  spends  about 
10%.  Enhancing  the  performance  of  the  moments  function  and  splitting/assembly  routines 
represents  the  largest  potential  gain  in  terms  of  increased  computational  efficiency 

The  execution  time  for  the  SC,  CLC,  and  LC  methods  is  fairly  constant.  The 
execution  speed  of  the  CEC  and  in  particular  the  EC  method  is  subject  to  variations  based 
on  several  factors.  First,  as  detailed  in  Appendix  A,  the  moment  function  routines  use 
either  recursion,  function  expansion  approximation,  or  a  combination  of  both  to  compute 
multi-argument  moment  functions.  The  expansions  execute  faster  than  recursion,  so  the 
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execution  time  of  a  multi-dimension  moment  function  will  vary  based  on  its  arguments. 
Second,  as  explained  in  Appendix  B,  the  root-solving  routines  can  use  two  separate 
formulations,  based  on  the  values  of  the  arguments,  to  determine  the  expansion 
coefficients.  The  single-argument  formulation  for  the  two-  and  three-dimensional  root 
solvers  is  more  efficient  than  the  multi-dimensional  form,  affecting  the  problem  run  time. 
Finally,  the  EC  method  uses  either  Equations  (88)  and  (90)  or  Equations  (108)  to 
compute  the  cell  flux  moments.  Using  the  balance  equations  is  more  efficient  than 
Equations  (88)  and  (90),  causing  variations  in  the  execution  speed  of  the  EC  method. 
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5.2.  Test  Problem  2  -  Cube  Problem  with  Isotropic  Incident  Current  and  Slab 
Symmetry 

Test  problem  2  consists  of  a  homogeneous  unit  cube  with  reflective  boundaries  at 
x  =  ±  0.5  cm  and  y  =  ±  0.5  era,  and  vacuum  boundaries  at  z  =  ±  0.5. cm  An  isotropic 
current,  Jz(x,y,  -0.5  cm)  =  1.0  cnT2  sec"1 ,  is  incident  on  the  z  =  -0.5  cm  face.  The  flux 
produced  by  the  isotropic  entering  current  is  a  function  of  the  streaming  angle,  Cin : 


2  J*(x,y,  -0.5  cm) 

-0  5cm)  =  — - - - - ,  (£„>0). 


(117) 


An  isotropic  incident  current  would  be  produced  by  an  isotropically  emitting  surface- 
deposited  source.  This  boundary  condition  is  used  since  it  is  more  stressing  than  a 
Lambertian  boundary  condition.  The  Lambertian  incident  current  consists  of  a  constant 
entering  angular  flux,  i.e.  i|rn  =  constant,  which  can  give  unexpectedly  good  results  that  are 
not  a  true  measure  of  the  quadrature  performance.  The  cube  material  has  a  scattering 
cross  section  of  4  cm'1  and  a  total  cross  section  of  32  cm'1 .  The  meshes  for  this  problem 
were  all  generated  using  IDEAS  and  contained  from  6  to  22028  cells.  The  average  linear 
cell  thickness  for  this  problem  is  1 8  mfp  on  the  coarsest  mesh  and  1 . 1  mfp  on  the  finest 
mesh.  The  maximum  linear  cell  dimension  is  54  mfp  for  the  coarsest  and  3  mfp  for  the 
finest  mesh.  The  TetSN  solutions  for  this  problem  were  generated  using  an  S8  level 
symmetric  quadrature  and  a  convergence  tolerance  of  lO*6. 

Because  the  problem  is  symmetric  about  two  axes,  we  can  compare  the  TetSN 
results  to  a  one-dimensional  semi-analytic  solution  with  equivalent  boundary  conditions 
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(Davison,  1958:  174-175).  The  semi-analytic  results  were  computed  in  Mathematica 
(Wolfram,  1992)  using  a  modified  three-dimensional  S8  level-symmetric  quadrature.  The 
modified  three-dimensional  quadrature  used  only  the  £  values  (direction  cosines  with 
respect  to  the  z-axis)  to  specify  the  angles.  The  weights  for  this  modified  quadrature  were 
just  the  sum  of  the  point  weights  of  all  the  angles  having  the  same  E,  value  (the  level 
weights). 

Plots  of  the  problem  average  flux  and  relative  error  of  the  average  flux  compared 
to  the  semi-analytic  solution  are  shown  in  Figure  28  and  Figure  29.  All  methods  are  able 
to  compute  the  average  scalar  flux  to  within  a  few  percent  even  on  the  coarsest  mesh,  with 
LC  and  in  particular,  EC,  doing  the  best  job.  Because  the  shield  is  optically  very  thick,  the 
flux  is  highest  near  the  entering  face  and  decreases  rapidly  in  the  +z  direction.  There  is  a 
thin  boundary  layer  near  the  entering  face,  where  the  entering  particles  are  scattered  back 
out  of  the  shield.  The  problem  average  flux  is  dominated  by  the  flux  in  this  region.  The 
constant  face  and  LC  methods  are  capable  of  adequately  representing  the  face  and  cell 
distributions  in  this  region  and  hence  do  a  reasonable  job  computing  the  problem  average 
flux.  Beyond  the  boundary  layer,  the  flux  decays  essentially  exponentially.  The  EC 
method  is  able  to  almost  exactly  model  the  flux  falloff,  and  computes  the  problem  average 
flux  to  within  about  1%  and  exiting  face  current  to  within  10%  using  only  six  cells.  The 
CEC  method,  because  it  is  better  able  to  model  the  scattering  source  distribution, 
performs  slightly  better  than  the  other  constant  face  methods  on  the  coarser  meshes,  but 
the  overall  performance  of  these  methods  is  poor. 
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Figure  28.  Average  Flux  for  Test  Problem  2. 
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The  most  stressing  calculation  is  the  exiting  face  (penetrating)  current.  The 
penetrating  current  and  relative  error  of  the  penetrating  current  compared  to  the  most 
converged  EC  solution  are  shown  in  Figure  30  and  Figure  31.  The  magnitude  of  the 
exiting  flux  computed  using  the  EC  method  is  essentially  constant  for  all  mesh 
refinements.  The  LC  solution  oscillates  between  positive  and  negative  values  as  the  mesh 
is  refined,  and  only  approaches  the  EC  result  on  the  finest  mesh.  Examining  the  relative 
error  plot,  we  see  that  the  LC  solution  converges  toward  the  correct  result,  but  even  the 
finest  mesh  LC  result  is  not  as  accurate  as  the  six  tetrahedron  EC  solution.  The  exiting 
current  predicted  by  the  constant  face  methods  differed  by  many  orders  of  magnitude  on 
the  coarser  meshes  and  is  still  off  by  two  orders  of  magnitude  on  the  finest  mesh. 


Figure  30.  Exiting  Face  Current  for  Test  Problem  2. 
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On  the  coarsest  mesh,  six  tetrahedra,  the  relative  errors  in  the  EC  solutions  for  the 
problem  average  flux  and  exiting  current  are  less  than  LC  solution  on  the  finest  mesh 
(22028  tetrahedron).  So  even  though  the  LC  solution  costs  1/8  as  much  per  cell  as  the 
EC  solution,  the  LC  method  requires  3375  times  as  many  cells  to  approach  the  same 
accuracy.  For  this  type  of  problem,  the  EC  method  is  clearly  the  most  computationally 
efficient. 
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Figure  3 1 .  Relative  Error  in  Exiting  Face  Current  for  Test  Problem  2. 


On  the  coarse  meshes,  LC  gives  unphysical  negative  exiting  currents.  The  fact  that 
the  CLC  method  remains  positive  indicates  that  it  is  the  face  distribution  that  causes  the 
LC  method  to  generate  negative  exiting  currents.  This  demonstrates  the  idea  that  CLC  is 
essentially  a  negative  face  flux  fix  up  for  LC.  Nevertheless,  the  LC  solutions  are  more 
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accurate  than  the  CLC  solutions,  which  raises  the  issue  of  positivity  versus  accuracy.  The 
LC  quadrature  produces  local  negative  values,  but  the  integral  quantities,  problem  average 
flux  and  exiting  face  current,  are  more  accurate  than  the  strictly  positive  (SC,  CEC) 
constant  face  method  solutions.  Further,  exiting  face  flux  and  current  values  produced  by 
the  constant  face  methods  that  differ  by  several  orders  of  magnitude  from  the  correct 
value  are  no  more  useful  than  the  negative  values  produced  by  the  LC  method  on  those 
meshes.  In  fact,  the  negative  LC  values  may  be  less  problematic  because  they  are 
obviously  nonsense,  whereas  the  constant  face  results  could  be  interpreted  as  being 
significant.  Strict  positivity  in  and  of  itself  is  not  necessarily  a  great  advantage  in  the 
absence  of  other  desirable  features.  Any  advantages  gained  by  the  strictly  positive 
behavior  of  the  constant  face  methods  in  this  problem  are  offset  by  the  increased  accuracy 
of  the  first-moment  conserving  LC  quadrature,  while  SC  and  LC  are  overshadowed  by 
EC,  which  is  both  accurate  and  strictly  positive. 

Table  3  compares  the  benchmark  and  the  finest  mesh  (22028  cell)  TetSN  EC 
solution  values  for  several  integral  quantities.  The  agreement  between  the  solutions  is 
very  good;  all  listed  quantities  agree  to  at  least  three  digits.  Although  the  differences  are 
small,  they  are  significant  since  the  TetSN  solutions  were  converged  to  six  digits.  The 
errors  in  the  TetSN  solution  are  most  likely  caused  by  numerical  diffusion.  Carrying  first 
moments  on  faces,  which  the  first  moment  methods  do,  reduces  numerical  diffusion,  but 
does  not  eliminate  it.  The  entering  current  computed  by  TetSN  is  slightly  higher  than  the 
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benchmark,  meaning  fewer  particles  were  back  scattered  at  the  entering  face.  A  small 
fraction  of  the  particles  are  numerically  diffusing  into  the  problem  interior.  Because  fewer 
particles  are  scattered  back  at  the  face,  the  TetSN  entering  face  flux  is  lower,  and  problem 
average  flux  and  exiting  face  flux  and  current  are  higher,  than  the  benchmark  solution. 


Table  3.  Comparison  of  One-Dimensional  Semi- Analytic  and  TETSN  Results  for 
Integral  Quantities  of  Test  Problem  2. 


Quantity 

ID  Semi- Analytic 
Result 

TETSN  Result 

Relative 

Error 

<i>A 

0.034523 

0.034545 

6.37x1 0^ 

4>(0) 

2.8513 

2.84097 

3.62x1 0‘3 

4>0) 

5.8095xl0"16 

5.8133xl0"16 

6.54X10"4 

ret  (o) 

0.96662 

0.96725 

6.52X10-4 

Fet  (1) 

5.2159xl0‘16 

5.2181xl0'16 

4.22x1  O'4 

Note  that  since  the  flux  is  decreased  by  16  orders  of  magnitude  in  traversing  the 
shield,  only  one  in  1016  particles  started  in  a  Monte  Carlo  solution  would  reach  the  exiting 
side  and  contribute  to  a  tally.  However,  because  of  the  symmetries  involved,  certain 
variance  reduction  techniques  could  be  employed  that  would  generate  an  efficient  and 
accurate  Monte  Carlo  solution  for  the  test  problem  2  exiting  flux  with  a  modest  number  of 
histories.  Difficulties  arise  when  spatial  information  about  fluxes  and  currents  are 
required,  especially  if  there  are  no  symmetries  to  exploit.  If  spatial  distributions  are 
needed  then  separate  estimators  (tallies)  are  needed  at  each  point  information  is  desired. 
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Each  new  estimator  added  to  a  computation  increases  the  computational  effort.  Further, 
the  variance  reduction  techniques  appropriate  for  one  estimator  are  not  necessarily  helpful, 
and  in  fact  may  be  detrimental,  in  reducing  the  variance  of  other  estimators  in  the  problem. 
For  more  complicated  structures  in  which  spatial  distributions  are  required,  the 
computation  costs  may  rapidly  become  prohibitive. 
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5.3.  Test  Problem  3  -  Cylinder-in-a-Cube,  Isotropic  Cylinder  Source 

Test  problem  3  consists  of  a  2x2x2  cm  cube  with  total  cross  section  of  4  cm'1  and 
scattering  cross  section  of  2  cm'1 .  Embedded  in  the  center  of  the  cube  is  a  0.5  cm  radius, 
1  cm  long  cylinder.  The  cylinder  material  has  a  total  cross  section  of  1  cm'1  and  a 
scattering  cross  section  of  0.5  cm'1 .  In  addition,  the  cylinder  material  contains  an 
embedded  isotropic  source  of  strength  1  cm0  sec'1 .  A  cutaway  view  of  the  2511 
tetrahedron  mesh  for  test  problem  3  is  shown  in  Figure  32.  The  TetSN  solutions 
generated  in  this  section  use  an  S16  level  symmetric  quadrature  and  a  convergence 
tolerance  of  10"6. 


Figure  32.  Cut  Away  View  of  the  Test  Problem  3 
25 1 1  Cell  Mesh. 
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Interpreting  the  convergence  for  this  problem  is  not  as  simple  as  in  the  previous 
cases.  For  any  problem  involving  curved  surfaces,  such  as  the  cylinder  in  this  problem, 
refining  the  mesh  also  changes  the  region  volumes.  When  IDEAS  meshes  a  curved 
surface,  it  picks  points  on  the  surface  of,  and  interior  to,  the  volume,  and  creates  cells  by 
connecting  those  points.  The  result  of  this  is  that  volumes  of  regions  with  concave 
boundaries  are  overestimated,  and  volumes  of  regions  with  convex  boundaries  are 
underestimated,  but  all  volumes  approach  the  correct  volume  as  the  mesh  is  refined.  It  is 
possible  to  compensate  for  this  by  repositioning  the  surface  ,  but  without  prior  knowledge 
of  exactly  how  IDEAS  is  going  to  generate  cells,  there  will  always  be  some  error.  Table  4 
lists  the  cylinder  volume  and  error  relative  to  the  exact  volume  for  all  mesh  refinements. 


Table  4.  Volume  of  Mesh  Representation  of  Cylinder  (Source  Region)  for 
Test  Problem  3  and  Error  Relative  to  the  Exact  Volume  (0.78539  cnr1). 


# 

Cylinder 

Cells 

# 

Cube 

Cells 

Total  Cells 

Cylinder 

Volume 

%  Error 

55 

284 

339 

0.69207 

11.88 

174 

731 

905 

0.73240 

6.74 

366 

1432 

1798 

0.75689 

3.63 

593 

1918 

2511 

0.76412 

2.71 

1272 

5497 

6769 

0.77238 

1.66 

1844 

9597 

11441 

0.77611 

1.18 
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The  problem,  cube  material,  and  cylinder  material  average  fluxes  are  shown  in 
Figure  33  through  Figure  35,  along  with  the  MCNP  results.  The  MCNP  uncertainty  for 
each  of  the  average  fluxes  is  0.02%.  Several  trends  are  apparent.  First,  in  each  plot,  the 
difference  between  the  first  order  and  constant  face  solutions  is  obvious.  The  constant 
face  methods,  as  they  have  in  every  test  problem,  converge  more  slowly  to  the  solution. 
However,  on  the  finest  mesh  the  constant  face  relative  error  with  respect  to  the  MCNP 
result  is  3%,  only  slightly  larger  than  the  1.3%  first  moment  method  error. 


Figure  33.  Problem  Average  Flux  for  Test  Problem  3. 
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Figure  34.  Cube  Material  Average  Flux  for  Test  Problem  3. 


Figure  35.  Cylinder  Material  Average  Flux  for  Test  Problem  3. 
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The  difference  in  the  constant  face  and  first  moment  method  solutions  is  larger  in 
Figure  36  and  Figure  37,  plots  of  the  exiting  currents  on  faces  parallel  and  perpendicular 
to  the  axis  of  symmetry  of  the  source  cylinder.  The  MCNP  uncertainty  for  each  exiting 
current  is  0.2%.  In  this  case,  the  EC  method  does  a  slightly  better  job  computing  the 
exiting  current  on  coarser  meshes,  but  both  LC  and  EC  reach  the  same  solution  on  the 
finest  mesh.  Again,  the  constant  face  solutions  converge  much  more  slowly.  The  finest 


Figure  36.  Exiting  Current  at  z  =  -1  cm  for  Test  Problem  3. 


mesh  constant  face  solution  for  the  exiting  face  flux  has  about  the  same  error  relative  to 
the  MCNP  solution,  15%,  as  the  coarsest  mesh  EC  result.  To  make  any  comparison  of 
the  TetSN  and  MCNP  solutions  at  any  one  refinement  is  difficult  because  of  the  source 
volume  error  in  the  TetSN  mesh.  However,  as  the  mesh  is  refined,  all  solutions  approach 
the  MCNP  solution,  indicating  that  the  code  is  operating  correctly. 
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Figure  37.  Exiting  Current  at  y  =  -1  cm  for  Test  Problem  3. 


The  cylindrical  source  drives  all  the  fluxes  and  currents  in  the  problem.  Part  of  the 
error  in  the  TetSN  solution  at  each  mesh  point  is  due  to  the  error  in  the  source  volume.  In 
essence,  because  the  source  volume  changes,  each  mesh  refinement  generates  a  slightly 
different  problem.  This  masks  the  true  behavior  of  the  quadratures.  Corrected  values  for 
the  fluxes  and  currents  can  be  computed  by  multiplying  the  TetSN  results  by  the  ratio  of 
the  correct  cylinder  volume  to  the  meshed  cylinder  volume.  Plots  of  both  the  corrected 
and  uncorrected  values  for  the  problem  average  flux  and  exiting  current  are  shown  in 
Figure  38  and  Figure  39.  The  relative  errors  in  the  corrected  solutions  are  much  lower 
than  in  the  uncorrected  case.  The  worst  case  error  in  the  first  order  solution  for  the 
average  flux  drops  from  about  13%  to  3%,  and  the  worst  case  error  in  the  exiting  current 
drops  from  23%  to  about  12%. 
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It  is  not  clear  that  this  volume  correction  approach  works  in  general.  If  we  were 
to  plot  corrected  values  of  the  exiting  flux  computed  with  the  constant  face  methods, 
shown  in  Figure  37  and  Figure  36,  the  relative  errors  would  actually  increase.  A  more 
consistent  approach  may  be  to  position  the  problem  surfaces  so  that  the  meshed  volumes 
better  approximate  the  actual  region  volumes.  However ,  as  already  pointed  out,  this 
requires  some  advance  knowledge  of  how  IDEAS  will  generate  the  cells,  which  may  not 
always  be  obvious.  Further,  this  approach  may  prove  to  be  very  tedious  because  the 
surfaces  have  to  be  repositioned  before  generating  the  cells  for  each  refinement.  A  better 
solution  to  volume  correcting  the  results  may  be  to  devise  a  mesh  generator  that 
automatically  conserves  volumes. 


Figure  38.  Problem  Average  Flux  for  Test  Problem  3.  *  Indicates  Results  are 
Corrected  for  the  Source  Volume. 
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Figure  39.  Exiting  Current  at  z  =  -1  cm  for  Test  Problem  3.  *  Indicates  Results 
are  Corrected  for  Source  Volume. 
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5.4.  Test  Problem  4  -  Cylinder-in-a-Cube,  Isotropic  Incident  Current 

Test  problem  4  consists  of  an  infinite  slab  cell  with  a  periodic  array  of  embedded, 
cylindrical  absorbers.  The  calculation  is  performed  on  a  unit  cell  of  the  array,  giving  the 
same  cylinder  in  a  cube  structure  used  in  test  problem  3  but  with  different  boundary 
conditions.  The  isotropic  source  in  the  cylinder  is  replaced  by  an  isotropic  current  incident 
at  z=-l  cm.  The  cube  material  has  a  scattering  cross  section  of  7  cm'1  and  a  total  cross 
section  of  10  cm'1 .  The  cylinder  is  a  strongly  absorbing  material  with  scattering  cross 
section  of  2  cm'1  and  total  cross  section  of  20  cm'1 .  The  symmetry  boundaries  exist  at 
x  =  ±1  cm  and  y  =  ±1  cm ,  vacuum  boundaries  exist  at  z  =  ±1  cm. 

This  problem  is  a  modified  version  of  the  deep  shield  problem  evaluated  in  test 
problem  2.  Since  we  already  know  that  the  constant  face  methods  perform  poorly  on  this 
type  of  problem,  we  concentrate  on  the  first  moment  methods.  The  problem  average  flux 
computed  using  the  LC  and  EC  spatial  quadratures  is  shown  in  Figure  40.  The  average 
flux  in  the  cube  and  cylinder  materials  is  shown  in  Figure  41  and  Figure  42.  In  all  cases, 
the  EC  result  is  essentially  converged  on  the  coarsest  mesh.  The  LC  quadrature  is  able  to 
achieve  results  accurate  to  within  6%  of  the  converged  result  for  the  problem  and  cube 
material  average  fluxes  on  fairly  coarse  meshes.  This  is  again  due  to  the  average  flux 
being  dominated  by  the  boundary  layer  near  the  entering  face.  The  flux  away  from  the 
boundary  layer  contributes  only  a  small  percent  to  the  average  flux.  Since  the  EC  method 
better  approximates  the  actual  distribution  deep  within  the  shield  even  on  coarse  meshes,  it 
is  able  to  correctly  account  for  the  deep  shield  contribution  for  all  mesh  refinements. 
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Figure  4 1 .  Problem  Average  Flux  for  Test  Problem  4. 


Figure  40.  Cube  Material  Average  Flux  for  Test  Problem  4. 
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The  LC  method  does  not  perform  as  well  in  computing  the  cylinder  material 
average  flux  until  the  mesh  is  fairly  refined.  While  the  coarsest  mesh  (339  cells)  LC  results 
for  the  problem  and  cube  average  flux  differ  from  the  fully  converged  EC  result  by  about 
6%,  the  method  required  11441  cells  to  reach  the  same  level  of  accuracy  for  the  cylinder 
average  flux.  This  is  again  because  the  particles  contributing  to  the  cylinder  flux  have  to 
traverse  part  of  the  cube  material,  where  they  are  strongly  attenuated.  Upon  reaching  the 
cylinder,  they  are  attenuated  even  more  strongly.  There  is  a  boundary  layer  at  the  leading 
edge  of  the  cylinder  material,  but  it  is  not  as  large,  nor  as  significant  as  the  entering  face 
boundary  layer.  For  the  LC  method  to  compute  the  cylinder  flux  correctly,  the  mesh  must 
be  fine  enough  to  allow  the  method  to  accurately  model  the  flux  fall  off  deep  in  the  shield. 


Figure  42.  Cylinder  Material  Average  Flux  for  Test  Problem  4. 
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The  most  stressing  calculation  is  the  exiting  face  current.  The  average  exiting  face 
current  computed  by  the  quadratures  at  each  refinement  is  plotted  in  Figure  43.  Again, 
the  EC  method  is  essentially  converged  on  the  coarsest  mesh,  while  the  LC  method 
overestimates  the  current  by  several  orders  of  magnitude.  The  LC  result  still  shows 
significant  error,  12%,  on  the  finest  mesh.  The  main  difference  between  this  problem 
and  test  problem  2  is  the  presence  of  the  absorber.  In  test  problem  2,  the  exiting  flux 
distribution  was  constant.  In  this  case,  the  exiting  flux  has  a  nonuniform  spatial 
distribution  due  to  the  presence  of  the  absorber.  Computing  the  exiting  flux  distribution  is 
especially  difficult  in  this  case  because  the  flux  is  already  highly  attenuated  by  the  shield 


Figure  43.  Exiting  Face  Current  for  Test  Problem  4. 
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(cube)  material.  Without  the  cylindrical  absorber,  the  flux  would  decrease  many  orders  of 
magnitude  in  traversing  the  shield.  The  addition  of  the  cylindrical  absorber  makes  a  small, 
but  significant,  change  in  the  exiting  flux  distribution. 

The  exiting  flux  distributions  computed  using  the  EC  and  LC  methods  for  each 
refinement  are  plotted  in  Figure  44  and  Figure  45.  The  EC  quadrature  senses  the  presence 
of  the  absorber,  and  has  essentially  the  correct  distribution  on  the  coarsest  mesh.  On  the 
finest  mesh,  the  EC  method  smoothly  models  the  distribution,  and  the  slope  of  the 
distribution  on  the  boundary  (since  the  problem  contains  an  infinite  array  of  embedded 
cylinders,  the  derivative  of  the  face  flux  at  x=±l  and  y=±l  must  be  zero).  The  LC  method 
is  only  able  to  resolve  the  flux  dip  on  the  finest  mesh,  and  is  never  able  to  get  the  slope  on 
the  face  edges. 

Figure  44  and  Figure  45  show  the  flux  values  plotted  using  the  mesh  nodes 
defining  each  face,  and  the  flux  distribution  coefficients  for  the  assumed  form  of  the 
distribution.  For  the  linear  distribution,  this  shows  all  the  available  detail.  For  the  EC 
solution,  the  solution  on  each  face  is  exponential,  and  there  is  detail  internal  to  each  face 
as  shown  in  Figure  46.  This  shows  that  EC  does  an  even  better  job  of  modeling  the  face 
flux  on  the  coarse  mesh  than  Figure  44  would  indicate. 
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Figure  44.  Test  Problem  4  Exiting  Face  Flux  Distribution  for  the  EC 
Method  (vertical-axis  is  flux  *  2x1 07) 


Figure  45.  Test  Problem  4  Exiting  Face  Flux  Distribution  for  the  LC 
Method  (vertical-axis  is  flux  *  2x1 07) 


131 


Figure  46.  Test  Case  4  Exiting  Flux  Distribution  Computed  Using  EC 
on  Two  Coarsest  Meshes  Using  Nodes  Values  (Left  Side)  and  Using 
Internal  Data  Points  (Right  Side),  (vertical-axis  is  flux*2xl07) 

A  quantitative  analysis  of  the  convergence  of  the  exiting  face  distribution  was  done 
using  the  error  norm  analysis  described  in  Section  5.1.1.  The  LOT  and  Lj  norms  of  the 
errors  in  the  exiting  face  distribution  are  shown  in  Figure  47  and  Figure  47.  While  the 
maximum  error  in  the  LC  solution  decreases  rapidly  with  decreasing  cell  size,  the  EC  error 
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for  each  refinement  is  still  significantly  less.  The  average  error  in  the  coarsest  mesh  EC 
solution  is  less  than  the  finest  mesh  LC  solution.  In  this  case  the  LC  method  requires 
more  than  37  times  as  much  computational  effort  to  achieve  the  same  accuracy  in  the  L} 
norm  as  the  EC  method.  This  demonstrates  the  advantages  of  the  EC  method  in  deep 
shielding  problems  for  pointwise  as  well  as  integral  measures. 


E 


-o-  EC  LC 


Figure  47.  Maximum  Pointwise  Relative  Error  of  the  Exiting  Face  Flux 
Distribution  for  Test  Problem  4. 


Again,  as  in  test  problem  2,  a  Monte  Carlo  method  would  have  a  very  difficult  time 
obtaining  results  for  the  average  exiting  current  because  of  the  probability  of  the  particles 
traversing  the  shield.  Computing  the  exiting  current  distribution  is  even  more  stressing 
because  the  Monte  Carlo  method  would  then  have  to  bin  the  surface  into  many  sub- 
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surfaces.  Then  each  sub-surface  would  require  as  many  particles  to  compute  its  average 
current  as  the  entire  face  required  to  compute  the  average. 


-a-  EC  -a-  LC 


Figure  48.  Average  Pointwise  Relative  Error  of  the  Exiting  Face  Flux 
Distribution  for  Test  Problem  4. 
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VI.  Summary/Recommendations 


Characteristic  spatial  quadratures  for  performing  discrete  ordinates  calculations 
using  arbitrary  tetrahedral  meshes  were  derived,  implemented,  and  tested.  Tetrahedral 
meshes  offer  several  advantages  over  traditional  boxoid  meshes.  First,  inclined  or  curved 
surfaces  can  be  accurately  modeled  using  fewer  cells.  More  importantly,  meshes  can  be 
refined  in  localized  areas  of  the  problem  without  drastically  affecting  the  cell  sizes  in  the 
rest  of  the  problem,  as  is  required  for  rectangular  parallelpiped  cells. 

The  algorithm  used  a  split  cell  approach  that  reduced  the  BTE  to  a  single 
dimension,  simplifying  the  characteristic  solution  and  the  resulting  moments  equations. 
Implementation  of  the  algorithm  required  that  the  methods  be  derived  for  the  special  case 
tetrahedron  having  a  single  inflow  and  a  single  outflow  face.  This  included  developing 
routines  that  could  numerically  evaluate  exponential  moments  functions  with  four 
arguments.  Methodology  for  defining  and  passing  moments  of  the  face  flux  distribution 
were  also  developed.  For  the  EC  method,  two-  and  three-dimensional  root  solvers  were 
developed  using  both  single  and  multiple  argument  formulations  to  increase  efficiency 
while  assuring  stability.  Equations  to  provide  an  initial  guess  which  guaranteed 
convergence  of  Newton' s  method  in  fewer  than  four  (generally  fewer  than  three) 
iterations  for  the  face  flux  and  source  distribution  systems  were  also  developed. 

Monte  Carlo  methods  have  become  the  preferred  method  for  doing  general- 
geometry  three-dimensional  neutral  particle  transport  due  both  to  the  inherent  benefits  of 
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the  Monte  Carlo  methods  and  the  lack  of  an  efficient  general-geometry  3D  discrete 
ordinates  code.  Aside  from  some  special  purpose  codes  designed  for  specific  applications, 
3D  discrete  ordinates  codes  are  limited  to  low  order  spatial  quadratures  and  rectangular 
parallelpiped  cells.  Although  there  are  many  problems  for  which  Monte  Carlo  methods 
are  an  appropriate  choice,  there  are  many  situations,  particularly  when  the  media  are 
optically  thick  or  when  spatial  distributions  are  required,  when  discrete  ordinates  may  be 
more  computationally  efficient.  For  what  ever  the  reason,  the  development  of  improved 
spatial  quadratures  has  not  been  translated  into  the  development  of  improved  operational 
codes,  particularly  in  three  dimensions.  The  TetSN  code  is  intended  to  be  a  first  step 
toward  the  development  of  a  general-geometry,  three-dimensional  code  package  that  will 
provide  a  needed  computational  tool  to  the  transport  community. 

6.1  Summary  of  Results 

For  test  problems  1,  2,  and  3,  the  characteristic  quadrature  solutions  converged  to 
the  independent  benchmark  result,  as  both  the  spatial  and  angular  variables  were  refined. 
Both  the  first  moment  methods,  LC  and  EC,  exhibited  third  order  convergence,  consistent 
with  the  convergence  rates  observed  for  the  methods  in  one  and  two  dimensional 
geometries.  The  constant  face  methods  were  all  approximately  first  order  convergent, 
again  consistent  with  the  SC  method's  performance  in  other  geometries.  The  relatively 
slow  convergence  of  the  constant  face  methods  on  fine  meshes,  and  the  resulting  low 
accuracy  and  severe  numerical  diffusion  on  coarse  meshes,  severely  limits  their  use  in 
practice.  Based  on  the  cases  examined,  there  appears  to  be  no  benefit  whatsoever  of  using 
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a  constant  face  flux  approximation  with  a  higher  order  source  distribution  approximation. 
The  SC  quadrature  performed  as  well  as  the  CLC  and  CEC  methods  on  each  of  the  test 
problems. 

All  of  the  quadratures  proved  to  be  relatively  insensitive  to  variations  in  the  mesh 
structure.  The  results  show  that  relative  to  the  precision  of  the  computed  results,  the 
fluctuations  in  solution  caused  by  variations  in  the  mesh  structure  were  essentially 
insignificant.  This  is  a  very  important  quality.  The  meshes  generated  in  Section  5.1.2, 
contained  some  very  poorly  shaped  tetrahedra.  In  some  situations,  it  may  be  necessary  to 
model  an  oddly  shaped  region  using  tetrahedra  with  poor  aspect  ratios.  For  these 
problems  it  is  essential  that  the  computation  method  be  insensitive  to  variations  in  the 
shape  of  the  cells. 

In  terms  of  computation  time  per  phase  space  cell,  the  EC  method  was  the  most 
expensive,  requiring  13  times  the  effort  of  SC  and  eight  times  the  effort  of  LC.  However, 
the  EC  method  excelled  on  problems  with  optically  thick  cells.  For  test  problem  2,  the  EC 
method  required  1/675  of  the  effort  of  the  LC  method  to  achieve  the  same  accuracy.  In 
test  problem  4,  where  the  cells  were  not  as  thick,  the  EC  method  required  1/37  of  the 
computational  effort  of  LC  to  compute  the  integral  quantities,  while  the  LC  method  never 
reached  the  coarse  mesh  EC  accuracy  in  computing  the  exiting  flux  distribution. 

On  problems  where  the  cell  optical  thickness  was  in  the  lmfp  range,  the  LC 
method  was  the  best  performer.  In  test  problems  1  and  3,  the  EC  and  LC  methods  gave 
about  the  same  results  at  each  level  of  mesh  refinement,  but  LC  required  only  1/8  of  the 
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EC  effort  to  achieve  the  same  accuracy.  However,  the  fact  that  the  EC  method  provides 
accurate  results  on  both  thin  and  thick  cells  is  a  very  desirable  trait.  Since  in  general 
multi-group  problems,  the  optical  thickness  of  a  material  can  change  several  orders  of 
magnitude  based  on  the  particular  energy  of  the  particles,  it  is  important  to  have  a  method 
that  is  stable  regardless  of  the  cell  size.  In  this  respect  the  EC  method  was  again  the  best 
performer. 

The  advantages  of  discrete  ordinates  over  Monte  Carlo  solutions  for  optically  thick 
problems  were  also  considered.  For  Problems  2  and  4,  where  the  flux  dropped  many 
orders  of  magnitude  in  traversing  the  shield,  the  Monte  Carlo  method  would  have  required 
a  prohibitive  number  of  particles  to  compute  the  average  results  without  considerable 
variance  reduction  assistance.  To  compute  distributions,  the  Monte  Carlo  method 
requires  many  times  more  particles  than  are  required  to  compute  the  average  quantities, 
further  compounding  the  difficulty.  The  flux  and  current  distributions  are  naturally 
computed  as  part  of  the  discrete  ordinates  solution. 

Finally,  the  advantages  of  using  tetrahedral  meshes  were  observed.  Just  as  in  the 
triangle  case  for  two  dimensions,  tetrahedral  cells  are  better  able  to  model  curved  surfaces. 
This  is  evidenced  in  test  problems  3  and  4  where  the  tetrahedra  were  able  to  model  the 
volume  and  general  shape  of  the  central  cylindrical  region  with  a  relatively  small  number 
of  cells.  This  will  prove  very  beneficial  in  problems  with  many  curved  surfaces. 
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6.2  Recommendations  for  Future  Efforts 


In  order  for  the  code  to  be  useful  as  a  production  tool,  several  limitations  of  the 
current  implementation  must  be  addressed.  Most  important  is  execution  speed.  Running 
the  test  problems  required  anywhere  from  several  seconds  to  over  a  day  to  complete.  In  a 
production  environment,  this  is  unacceptable.  The  current  implementation  was  intended 
to  be  straightforward  and  understandable.  There  are  doubtless  may  opportunities  for 
optimization  including  coding  practices,  compiler  optimizations,  and  even  operating 
system  optimizations.  The  area  in  which  modest  investment  of  effort  would  pay  the  most 
dividends  is  in  increasing  the  efficiency  of  the  exponential  moments  functions  routines. 
Currently  the  EC  and  LC  methods  spend  about  60%  and  30%  respectively  of  their 
execution  time  in  the  moments  function  routines.  All  moment  function  evaluations 
eventually  reduce  to  computing  a  number  of  single-argument  functions.  The  single 
argument  routines  were  written  to  be  able  to  accurately  compute  functions  of  any  order. 
However,  as  can  be  seen  from  the  equations  in  Chapter  3,  the  overwhelming  majority  of 
the  functions  required  are  only  of  order  0  and  1 .  Currently,  the  routines  for  numerically 
evaluating  the  moments  functions  are  able  to  treat  all  orders.  By  writing  routines  to  deal 
only  with  zeroth  and  first  order  functions,  it  may  be  possible  to  gain  considerable 
computational  efficiency. 

For  the  EC  method,  it  may  be  possible  to  gain  some  computation  speed  by  defining 
a  new  set  of  functions.  As  shown  in  Chapter  3,  all  the  equations  for  the  EC  method  first 
moments  are  in  terms  of  the  ratio  of  a  first  order  to  a  zeroth  order  function  of  the  same 
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arguments.  As  shown  in  Appendix  B,  the  domain  of  these  functions  is  bounded  between 
zero  and  one.  Rather  than  trying  to  compute  the  moments  functions  in  the  numerator  and 
denominator,  it  may  be  more  efficient  to  treat  the  ratio  itself  as  a  function,  and  develop 
routines  to  numerically  evaluate  these  new  functions. 

Increasing  the  efficiency  of  the  root-solvers  is  also  necessary.  Currently,  the  time 
spent  in  the  root-solvers  is  only  a  fraction  of  the  total  EC  execution  time,  but  addition  of 
anisotropic  scatter  makes  the  rootsolver  dominant.  The  code  currently  assumes  only 
isotropic  scatter.  For  many  problems  of  interest  this  assumption  is  too  limiting.  Adding 
the  capability  to  handle  anisotropic  scatter  is  straightforward,  but  adds  considerable 
computational  effort.  Currently,  the  source  distribution  coefficients  only  have  to  be 
computed  once  per  iteration,  since  the  source  is  independent  of  angle.  Once  the 
assumption  of  isotropic  scatter  is  dropped,  the  source  becomes  a  function  of  angle,  and 
the  source  distribution  coefficients  have  to  be  computed  for  each  angle.  Referring  back  to 
Table  2,  Section  5.1.3,  for  an  S16  quadrature  (288  angles),  the  computation  of  the  source 
distribution  coefficients  for  the  EC  method  goes  from  being  almost  insignificant,  to 
dominating  the  computation. 

There  are  several  ways  that  the  root  solvers  might  be  optimized.  First,  the  initial 
guess  equations  can  be  refined.  Currently,  the  root-solvers  require  3.5  iterations  on 
average  to  converge.  Having  better  initial  guess  equations  would  reduce  the  required 
number  of  iterations.  It  may  even  be  possible  to  generate  a  set  of  equations  that  could 
generate  reasonable  approximations  to  the  coefficients  without  having  to  rootsolve. 
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Another  possibility  is  to  examine  alternate  numerical  techniques.  We  chose  Newton's 
method  because  it  is  quadratically  convergent,  and  because  we  had  the  capability  to 
analytically  evaluate  the  Jacobian.  There  are  methods  that  are  slightly  less  than 
quadratically  convergent,  that  use  a  secant  approximation  to  the  function  derivatives  that 
may  be  more  computationally  efficient. 

The  code  as  written  can  accommodate  multiple  energy  groups,  but  doing  so  adds 
considerable  computational  effort.  For  a  multi-group  problem,  the  solution  must  converge 
for  each  energy  group.  Therefore,  the  execution  times  listed  in  Table  2  have  to  be 
multiplied  by  the  number  of  energy  groups  in  the  problem  for  down  scatter  only,  and 
several  times  the  number  of  energy  groups  when  up  scatter  is  allowed.  This  underscores 
the  need  to  improve  the  efficiency  of  the  moment  function  and  root-solving  routines. 

As  the  efficiency  of  the  moment  function  and  root-solve  routines  increase,  the 
splitting  and  assembling  routines  will  become  a  larger  fraction  of  the  execution  time.  For 
the  non-exponential  methods,  splitting  and  assembling  already  represent  more  than  50%  of 
their  execution  time,  hence  improving  the  routines  could  considerably  improve  the 
performance  of  these  methods. 

Finally,  although  the  IDEAS  code  was  able  to  generate  the  tetrahedral  meshes 
used  in  the  test  problems  with  a  modest  level  of  manual  intervention,  more  complicated 
structures  will  be  quite  tedious  to  generate.  Regardless  of  the  shape  of  the  object,  the 
code  can  generate  a  tetrahedral  mesh  to  fill  it.  However,  when  the  geometry  includes  not 
only  complex  shapes,  but  several  regions  with  different  materials  as  well,  it  takes 
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considerable  manual  intervention  to  generate  the  mesh.  None  of  these  adjustments  are 
obvious  unless  the  operator  is  very  familiar  with  the  features  of  the  code.  Before  TetSN 
will  be  useful  as  a  production  tool,  a  more  user  friendly  approach  to  mesh  generation  must 
be  devised. 
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Appendix  A.  Exponential  Moment  Functions 

The  exponential  moment  functions  were  introduced  in  section  3.3.  This  Appendix 
presents  the  methods  used  to  numerically  evaluate  the  functions. 

A.l  One-Dimensional  Exponential  Moment  Functions 

Algorithms  for  numerically  evaluating  moments  functions  of  a  single  argument 
were  developed  by  Minor  (1993).  These  routines  are  the  basis  of  calculation  for  higher 
dimension  moment  functions. 

A.2  Two-Dimensional  Exponential  Moment  Functions 

Expressed  in  terms  of  recursion  on  number  of  arguments,  the  form  of  the  two- 
argument  moment  function  is 


M„(*i)-M„(s2) 

(*2  ~xi) 


(118) 


where  x2  *  xl.  In  general,  any  two  argument  moment  function  can  be  efficiently 
computed  using  this  definition  and  the  single  argument  routines  developed  by  Minor. 
However,  as  the  value  of  Xj  approaches  that  of  x2,  the  equation  becomes  numerically  ill 
conditioned.  When  this  occurs.  Equation  (118)  can  be  rewritten  as 


Mn(w  -  A w ,  w + Aw) 


Mn  (w  -  Aw)  -  MJw  +  Aw) 
2  Aw 


(119) 
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where w  =  (xt  +  x2  )/2  and  Aw  =  w  -  Xj .  In  this  form,  as  Xj -*x2,  w-0,  and  the  function  can 
be  approximated  using  a  series  expansion  about  Aw=0. 

In  view  of  the  orderless  property,  without  loss  of  generality,  any  first  partial 
derivative  of  an  exponential  moment  function  is  obtained  using 

q 

—Mn(Xl  =  -Mn(Xl  .  (120) 


For  functions  of  a  single  argument,  the  derivative  is  obtained  from  the  integral  definition, 
equation  (36): 

-^-Mn(x)  =Mn+1(X)  -Mn(x) .  (121) 


Higher  order  derivatives  of  a  single  argument  function  are  found  by  recursively  applying 
Equation  (121)  giving. 


dk 

dxk 


Mn(X)  = 


Ik- 1 


dx 


it- i 


d^_ 

dXk~l 


M„(x ) . 


(122) 


Using  Equation  (122)  the  expansion  of  Equation  (1 19)  to  O(Aw)6  is 

Mn(w  -  A w,w  +  Aw)  ~  Mn  (w)  -  Mn+ ,  (w) 

*  [M»  -  3 -  M,„W  ] 

+  [Mn{w)  -  5Mn+l(w)  +  10 Mn+2(w) 

^  10 1C,W  -  5M„4(».)  ]  4^  . 


(123) 
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The  expansion  was  evaluated  numerically  to  determine  the  required  order  of  expansion  for 
a  given  moment  function  order,  n,  and  variable  separation,  Aw.  The  criterion  used  was 
that  the  computed  value  contain  at  most  one  digit  of  accuracy  less  than  the  computation 
tolerance  used  to  compute  the  single-argument  functions.  Based  on  these  results,  the 
algorithm  for  computing  two-dimensional  moment  function  is  as  follows: 

If  (Aw  <  I0'3) 

Compute  moment  function  using  expansion  with  up  to  third  order  terms 

Else  If  (I  O'3  <  Aw  <  10'2) 

Compute  moment  function  using  expansion  with  up  to  fifth  order  terms 

Else 

Compute  moment  function  using  recursion  on  order  and  single-argument 

moments  functions  routines 

Endif 

The  expansions  and  algorithms  hold  without  regard  to  the  sign  of  the  individual 
arguments. 

If  the  variables  are  far  enough  apart  in  magnitude,  the  two-dimensional  moment 
function  can  be  computed  using  Equation  (118),  which  requires  that  M^x)  and  M^y) 
both  be  calculated  in  the  one-dimensional  routines.  If  the  variables  are  close  to  one 
another,  i.e.  Aw  is  small,  the  two-dimensional  function  is  evaluated  using  the  expansion. 
Numerically  evaluating  the  expansion  requires  computing  at  most  M^w)  to  Mj^w). 

This  turns  out  to  be  a  much  more  efficient  calculation  in  most  cases  because  the  one¬ 
dimensional  routines  only  need  to  compute  a  given  order,  then  use  recursion  to  get  the 
remaining  order  functions,  e.g.  compute  M^x)  then  use  recursion  to  get  M^w)  to 
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Computing  a  single  order  then  recursing  is  generally  faster  than  computing  two 
separate  functions,  so  the  moment  function  calculation  is  actually  faster  when  using  the 
expansion.  This  effect  is  even  more  pronounced  in  higher  dimension  functions. 

A.3  Three-Dimensional  Exponential  Moment  Functions 

Expressed  in  terms  of  recursion  on  number  of  arguments,  the  three-dimensional 
moment  functions  are 


K(X  l>x2’xi)  = 


M (Xj  ,x2  )  A-In  (x, ,  X3 ) 

(x3  ~X,) 


(124) 


where,  without  loss  of  generality,  the  variables  Xj ,  x2,  and  x3  are  arranged  in  ascending 
order.  In  general  any  three-dimensional  moment  function  can  be  computed  using  this 
definition  and  routines  capable  of  computing  two-dimensional  moment  functions. 
However,  numerical  ill  conditioning  occurs  when  the  value  of  x3  approaches  that  ofxt .  In 
those  situations,  an  alternate  form  of  Equation  (124),  given  by 


Mn(w  -  Aw,w + Ay,w  +  Aw) 


Mn(w-Aw,w+Ay)  -Mn(w  +  Ay,w  +  Aw) 
2  Aw 


(125) 


x,  +x3 

where  w  =  — - — ,  Aw  =  x3  -  w  =  w  -  x3 ,  and  Ay  =  x2  -  w.  Since  Ay  is  always  less  than 
Aw,  as  Aw  approaches  zero.  Equation  (125)  can  be  approximated  by  an  expansion  about 
Aw  =  0  and  Ay  =  0.  Expanding  about  Aw  =  0  and  Ay  =  0  and  only  keeping  terms  to 
order  five  in  Aw,  Ay,  or  products  of  Au  and  Ay  gives. 
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Mfw -A w,w  +A y,w  +Aw)  ~  -^Mn(w)  -  2Mn.I(w)  +Mn+2(w)j 


4v 

6 


-Mn(w)  +  3 Mn+l(w)  -  3 Mn+2(w)  +Mw+3(vv)]- 


24 


+  [  -KM  +  5K+1M  -  ioa^„+2(w)  + 10 Mn+3(w) 

Ay3  +Aw2A_y 


5M+4W+M+5(w)]. 


120 


+  [M„(w)  ~  6M+1(w)  +  15 Mn+2(w)  -  20Mn+3(w) 

Ay2  Aw2 


+  15 K+4(w)  -  6Mn+5(w )  +Mn+6(w)]- 


720 


[  -Mn(w)  +  7M„+1(w)  -  21M/i+2(w)  +  35M„+3(h>)  -  35*0) 

A^3A^2 


+  21M„+5(w)  -  7Mn+6(w)  +Mn+1{w)]- 


5040 


(126) 


As  for  the  two-dimensional  case,  the  expansions  were  evaluated  numerically  to  determine 
the  required  order  of  expansion  for  a  given  moment  function  of  order  n,  and  variable 
separation  Aw.  The  criteria  used  were  the  same  as  for  the  two-dimensional  case.  The 
algorithm  for  determining  which  method  to  use  in  computing  three-dimensional  moments 
functions  is: 

Sort  variables  in  ascending  order 
If  (Aw  <  10-3) 

Compute  moment  function  using  expansion  with  up  to  third  order  terms 
Else  If  (10 '3  <  Aw  <  1 0'2) 

Compute  moment  function  using  expansion  with  up  to  fifth  order  terms 
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Else 


Compute  moment  function  using  recursion  on  order 

Endif 

For  Aw  >  10'2,  recursion  on  order  gives  accurate  enough  results.  The  resulting  two- 
argument  functions  can  then  be  sent  to  the  two-dimensional  routines  for  evaluation. 
However,  if  the  separation  of  the  arguments  of  both  two-dimensional  functions  are  large 
enough  that  they  can  be  computed  by  recursion  to  first  kind  functions,  then  one  of  the  first 
kind  functions  ends  up  being  computed  twice  as  shown  in  Figure  49.  The  algorithm  for 
computing  three-argument  functions  tracks  the  variable  separation  at  each  step  and 
ensures  the  minimal  amount  of  computation. 


Figure  49  Details  of  the  number  of  single  moment  functions  that  have  to  be 
evaluated  to  compute  a  three  argument  function  by  recursion  on  number 
of  arguments. 
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With  this  in  mind,  the  algorithm  for  computing  three-argument  moment  functions  is 

Sort  Variables  in  ascending  order 

If  (x3-  x  f  <  I  O'2  then 

Compute  Mfxb  x2,  xf  using  expansion 

Else 

If  (x 2  -  Xj)  <  JO'2  and  (x3  -  xf  <  I  O'2  then 

Compute  Mn(xj,  xf  andMn(x2,  xf}  using  expansions 
Compute  Mf. xb  x2,  xf  using  recursion 
Else  If  (x 2  -  Xj)  <  I  O'2  and  (x3  -xf  >  I  O'2  then 
Compute  Mn(xb  xf  using  expansions 
Compute  Mfxf ,  andMfxf 
Compute  Mn(x2,  xf  using  recursion 
Compute  Mfxb  x2,  xf  using  recursion 
Else  If  (x 2  -  xf  >  10~2  and  (x3  -xf  <  10'2  then 
Compute  Mfxf,  Mfxf 
Compute  Mnfxj,  xf  using  recursion 
Compute  Mn(x2,  xf  using  expansion 
Compute  Mfxp  x2,  xf  using  recursion 
Else  if  (x2  -  xf  >  10~2  and  (x3  -xf  >  10'2  then 
Compute  Mfxf,  Mfxf,  and  Mfxf 
Mn(xi’  x2>  =  (Mfef  -Mfxf)/ (x2-x1} 

Mn(x2,  xf  =  (Mfxf  -Mfxf)/(x3  -  xf 
Compute  Mfxb  x2,  xf  using  recursion 

End  If 

End  If 

The  expansions  and  algorithm  are  not  dependent  on  the  sign  of  the  individual  arguments. 
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If  the  arguments  are  such  that  the  three-dimensional  function  can  be  computed 
using  an  expansion,  then  at  most  M^w)  to  M^^w)  will  have  to  be  computed.  Any  other 
situation  requires  that  at  least  two,  and  at  most  three,  separate  single  argument  functions 
be  evaluated.  Since  the  calculation  by  recursion  on  order  is  very  inexpensive,  evaluation 
of  the  function  using  the  expansion  tends  to  be  much  quicker  than  evaluation  by  recursion 
or  combination  of  recursion  and  expansion. 

A.4  Four-Dimensional  Exponential  Moment  Functions 

Expressed  in  terms  of  recursion  on  number  of  arguments,  the  four-dimensional 
exponential  moments  functions  are 


(Xj  ,x2,Xj)  A</j(x2,x3  ,x4) 

(*4~Xl) 


(127) 


where  the  variables  are  arranged  in  ascending  order.  In  general  any  moments  function 
having  four  arguments  can  be  computed  using  this  equation  and  routines  that  compute 
three-argument  moment  functions.  Numerical  ill-conditioning  occurs  as  the  value  of  Xj 
approaches  that  of  x4.  In  this  case  Equation  (127)  can  be  written  as, 

\in(yv  -  Aw,w  +  Av -Am,w  +  Av  +  Aw,w  +  Aw) 

Mn(yv  -  A>r,w  +  Av  -A u,w  +  Av  +A«)  ~Mn(w  +  Av  -A u,w  +  Av  +A«,w  +  Aw)  >0^9) 

— 


X4  +  Xj  X2  +  X3 

where  w  = - ,  v  = - ,  Aw  =  x4  -  w,  Av  =  v  -  w,  and  Au  =  x3  -  v.  Since  Au  and 
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Av  are  always  less  than  Aw,  as  Aw-O,  Equation  (128)  can  be  expanded  in  a  series  about 
Aw=0,  Av  =  0,  and  Au  =0.  Keeping  terms  to  order  five  in  Au,Av,Aw,  or  products  of 
Au,Av,  and  Aw, 


Mn(w  -  A w,w  +Av  -  A u,w  +  Av  +Aw,w  +  Aw) 

*  -  3K+l(W)  +  3Mn^(W) 


+  \-Mn(w)  +  4 Mn+X{w)  -  6Mn+2(w)  +  4Mn+3(w)  -  A/w+4(w)]-^ 

+  [Mn{w)  -  5 Mn+l(w)  +  \0Mn+2(w)  -  \0Mn+3(w) 

3  Av2  +  Aw2  +  At/2 


+  5Mn+4(w)  -Mn+5(w)]- 
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+  [  ~MJw)  +  6 Ma+1(w)  -  l5Mn+2(w)  +  20Mn+3(w)  -  15 Mn+4(w) 

2Av3  +  AvAw2  +2A«2Av 


6Mn+5(w)  -Mn+6(w)]- 


180 


[M„(w)-7Mn+l(w)+21Mn+2(w) 

~  35 Mn+3(w)  +  35 Mn+4(w)  -  21  Mn+5(w)  +  7M„+6(iv) 
3  Av2Aw2  +  Am  2  Aw  2  +  10Am2Av2 


5040 


+  [  -K(w)  +  *K+i(w)  - 2871 4n+2(w)  +  5 6M+3(w) 

-  70 M„+4(W)  +  56Mn+5(>v)  -  2SMn+6(w)  +  87V/„+» 

-  TV/  (V)  1  Av3Aw2  +  An2AvA>v2  +  5AM2Av3 
"+9  10080 


(129) 


For  Aw  >10'2,  then  functions  are  computed  using  recursion.  For  Aw  <10'2,  the  expansion 
using  up  to  fifth  order  terms  must  be  used,  and  for  Aw  <  10'3,  the  expansion  using  up  to 
third  order  terms  is  used  to  evaluate  the  functions. 
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If  the  arguments  are  such  that  the  four-argument  function  can  be  computed 
entirely  by  recursion  to  first  kind  functions,  care  must  be  taken  to  avoid  unnecessary 
duplication  of  effort.  Figure  50  shows  that  blindly  recursing  through  the  lower 
dimensional  routines  results  in  twice  as  many  first  kind  evaluations  as  needed.  The 
algorithm  for  computing  four-argument  functions  tracks  the  variable  separation  at  each 
step  and  ensures  the  minimal  amount  of  computation. 


Figure  50.  Number  of  single-argument  moment  functions  that  have  to  be  computed  to 
evaluate  a  four-argument  function  completely  by  recursion  on  number  of  arguments. 
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With  this  in  mind,  the  algorithm  for  computing  four-dimensional  moment  functions  is: 

Sort  variables  in  ascending  order 
If  (x 4-  Xj)  <  I  O'2  then 

Compute  Mn(X],  x2,  x3,  xf  using  expansion 
Else  If  (x3  -  x})  <  10~2  and  (x4  -  xf  <10~2  then 

Compute  Mn(Xj,  x2,  x^  andMn(x2,  x3,  xf  using  expansion 
Compute  Mn(xb  x2,  x3,  xf  using  recursion 
Else  If  (x3  -  x})  >  I  O'2  and  (x4  -  xf  <  I  O'2  then 
If  (x 2  -x})  <  10'2  then 

If  (x3  -  xf  <  I  O'2  then 

Compute  Mn(x2,  x3,  xf  using  expansion 
Compute  Mn(xh  xf  andMn(x2,  xf  using  expansion 
Compute  Mn(xj,  x2,  x^  using  recursion 
Compute  Mn(xj,  x2,  x3,  xf  using  recursion 
Else  If  (x3  -X2)  >  I  O'2  then 

Compute  Mn(x2,  x3,  xf  using  expansion 
Compute  Mn(xj,  xf  using  expansion 
Compute  Mfxf)  andMfxf 
Compute  Mn(x2,  xf  using  recursion 
Compute  Mn(xj,  x2,  xf  using  recursion 
Compute  Mn(xj,  x2,  x3,  xf  using  recursion 

End  If 

Else  If  (x 2  - x })  >10-2  then 

If  (x3  -  xf  <  10'2  then 

Compute  Mn(x2,  x3,  xf  using  expansion 
Compute  Mn(x2,  xf  using  expansion 
Compute  Mfx ;)  and  Mn(xf) 

Compute  Mn(xh  xf  using  recursion 
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Compute  Mn(xj,  x2,  x^  using  recursion 
Compute  MJxj,  x2,  x3,  xf  using  recursion 
Else  If  (x3  -xf)>  10'2  then 

Compute  Mn(x2,  x3,  xf  using  expansion 
Compute  M/xj),  Mn(xf},  and  Mfxf 
Compute  Mn(xj,  x-f  andMn(x2,  x^  using  recursion 
Compute  Mn(xj,  x2,  x^  using  recursion 
Compute  Mn(x},  x2,  x3,  xf  using  recursion 

Endlf 

Endlf 

Else  If  (x3  -  Xj)  <  I  O'2  and  (x4  -x2)>  I  O'2  then 

If  (x2  -  Xj)  <  10' 2  then 

If  (x4-x3)  <  10'2  then 

Compute  Mn(xj,  x2,  xf)  using  expansion 
Compute  Mn(x2,  xj  and  Mn(x3,  xj  using  expansion 
Compute  Mn(x2,  x3,  xf  using  recursion 
Compute  Mn(Xj,  x2,  x3,  xf  using  recursion 
Else  If  (x4  -x3)>  10'2  then 

Compute  Mn(xj,  x2,  x^  using  expansion 
Compute  Mn( x2,  x^  using  expansion 
Compute  Mfxf  and  MJxf 
Compute  MJx3,  xf  using  recursion 
Compute  Mn(x2,  x3,  xj  using  recursion 
Compute  Mn(xj,  x2,  x3,  xf)  using  recursion 

Endlf 

Else  If  (x2  -Xj)  >10'2  then 

If(x4-x3)  <  10'2  then 

Compute  Mn(xj,  x2,  x^  using  expansion 
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Compute  Mn( x3,  xf  using  expansion 
Compute  Mn(x 2)  andMJx $) 

Compute  Mn(x2,  x^  using  recursion 
Compute  Mn(x2,  x3,  xj  using  recursion 
Compute  Mn(xj,  x2,  x3,  x^)  using  recursion 
Else  If  (x4-x3)  >  JO'2  then 

Compute  Mn(xp  x2,  xf  using  expansion 
Compute  Mn(x 2),  MJxf,  andMJxf 
Compute  Mn(x2,  xf  andMn(x3,  xf)  using  recursion 
Compute  Mn(x2,  x3,  xf  using  recursion 
Compute  Mn(xj,  x2,  x3,  xf  using  recursion 

End  If 

End  If 

Else  If  (xn-xf  >I0'2  and  (xj-xf  >  10'2  then 
If  (x2  -X])  <  I  O'2  then 

If  (x3-xf  <  10'2  then 

If  (x4-xf)  <  I  O'2  then 

Compute  Mn(xp  xf,  Mn(x2,  xf,  Mn(x3,  xf  using 
expansion 

Compute  Mn(xj,  x2,  xf  and Mn(x2,  x3,  xf  using 
recursion 

Compute  Mn(xj,  x2,  x3,  xf  using  recursion 
Else  If  (xj-xf  >  I  O'2  then 

Compute  Mn(xh  xf  andMn(x2,  xf  using  expansion 
Compute  Mn(Xj,  x2,  x^  using  expansion 
Compute  MJxf  andMJxf 
Compute  Mn(x3,  xf  using  recursion 
Compute  Mn(x2,  x3,  xf  using  recursion 


155 


Compute  Mn(xj,  x2,  x3,  xf  using  recursion 

Endlf 

Else  If  (x3-x 2)  >  1  O'2  then 

If  (x^-xf  <  10'2  then 

Compute  Mn(xf  and  Mn(xf 
Compute  Mn(x2,  xf  using  recursion 
Compute  Mn(xj,  xf  andMn(x3,  xf  using 
expansions 

Compute  Mn(xj,  x2,  xf  andMn(x2,  x3,  xf  using 
recursion 

Compute  Mn(xj,  x2,  x3,  xf  using  recursion 
Else  If  (xj-xf  >  l  O'2  then 

Compute  Mn(xj,  xf  using  function 

Compute  Mn(xf,  Mn(xf,  Mn(xf 

Compute  Mn(x2,  xf  andMn(x3,  xf  using  recursion 

Compute  M/xj,  x2,  xf  cmdMn(x2,  x3,  xf  using 

recursion 

Compute  Mn(xh  x2,  x3,  xf  using  recursion 

Endlf 

Endlf 

Else  If  (x2  -xf  >  I  O'2  then 
If  (x3-x 2)  <  I  O'2  then 

If  (x^-xf  <  10'2  then 

Compute  Mn(x2,  xf},  Mn(x3,  xf}  using  expansion 
Compute  Mn(xj)  andMfxf 
Compute  Mn(xj,  xf  using  recursion 
Compute  Mn(x]t  x2,  xf  andMn(x2,  x3,  xf  using 
recursion 
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Compute  Mn(xj,  x2,  x3,  xf  using  recursion 
Else  If  (x^xf  >  10' 2  then 

Compute  Mn(x2,  xf)  using  expansion 
Compute  Mfxj),  Mn(xf},  Mn(xf,  Mn(xf 
Compute  Mn(xj,  xf  andMn(x3,  xf  using  recursion 
Compute  Mn(xp  x2,  x^  andMn(x2,  x3,  xf  using 
recursion 

Compute  Mn(xj,  x2,  x3,  xf  using  recursion 

End  If 

Else  If  (x3-x  t)  >  10' 2  then 
If  (xj-xf  <  1  O'2  then 

Compute  Mn(x3,  xf  using  expansion 
Compute  Mn(x}),  M/xf  MJxf,  Mfxf 
Compute  MJx {,  xf  andMn(x2,  xf  using  recursion 
Compute  Mn(xj,  x2,  x^  andMn(x2,  x3,  xf  using 
recursion 

Compute  Mn(Xj,  x2,  x3,  xf  using  recursion 
Else  If  (x4-xf)  >  JO'2  then 

Compute  Mn(x}),  MJxf,  Mfxf,  MJxf 
Compute  Mn(xlt  xf,  Mn(x2,  xj,  Mn(x3,  xf  using 
recursion 

Compute  MJxj,  x2,  xf)  andMn(x2,  x3,  xf  using 
recursion 

Compute  Mn(xj,  x2,  x3,  x4)  using  recursion 

End  If 

End  If 

Endlf 

End  If 
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Appendix  B.  Root  Solving 


B.l.  Newton's  Method 

The  equation 

=  \id  -  J  l  ( *oid)  F (\id)  »  (130) 


describes  the  Newton  iteration  for  solving  a  nonlinear  system  of  equations  of  the  form 


F(!)  =[/;(!)  ,/2(l),..  .,/„(!)]  =0, 

where 


(131) 


(132) 


5F; 


The  entries  of  the  Jacobian  matrix,  J,  are  given  as  J  =  — 

U  dxj 

The  method  starts  with  some  initial  guess  to  the  vector  x.  It  then  evaluates  F(x) 
and  J_1(x)  and  uses  these  to  compute  a  correction  to  the  initial  guess.  This  process  is 
repeated  until  either  the  function  values  are  zero  to  a  given  tolerance,  or  the  change  in  the 
x  values  from  the  last  iteration  is  less  than  the  prescribed  tolerance  (Press,  et.  al.  1992: 
372-375). 


The  method  displays  quadratic  convergence,  given  an  appropriate  starting  guess. 
If  the  initial  guess  is  poor,  the  method  may  take  many  iterations  to  converge,  or  may  not 
converge  at  all.  The  initial  guess  equations  developed  for  the  two-  and  three-dimensional 
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root  solvers  generate  first  guesses  that  generally  lead  to  convergence  to  the  desired 
tolerance  in  three  to  four  iterations. 

B.2.  Calculation  of  the  Input  Flux  Distribution  Coefficients  -  2D  System 

B.2.1.  Derivation  of  2D  System  Using  Multi-Dimensional  Moment  Functions 

As  explained  in  Section  3.4,  the  coefficients  of  the  exponential  distribution  must  be 
generated  from  knowledge  of  the  first  two  spatial  moments  of  the  distribution.  The 
moments  of  the  entering  flux  distribution  are  given  by 


c  =  2  aex»M0(Xin,YJ 

=  €  [Ap«»,4)] 

=  <  [y*p(-4,4)  ♦  ip«„-4,-4)]. 

where 

4 


(134) 


(134) 


Given  the  equations  for  the  flux  moments,  define 


_  < 

y  h c  KWM 


(136) 


and 
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i  in 


(136) 


y9y  M0(Xin-Yin,-Yiny 


Further,  define 


P^  =  Px-YP>. 


(137) 


Equations  (136)  and  (137)  define  a  system  of  two  equations  in  the  two  unknowns  and 
Y^.  Since  the  system  is  non-linear,  the  values  of  Xin  and  Yin  must  be  found  using 
Newton's  Method  on  the  system 


fiX.  ^(^4),n 

/i(  **  J  M0(Xin,YJ 


(138) 


and 


/,(«„)  =  P. 


xy 


(139) 


B.2.2.  Note  on  Notation 

The  routines  for  doing  the  two-dimensional  root  solve  were  originally  written  for 
the  triangle  code,  and  the  notation  was  slightly  different.  As  shown  in  Figure  5 1,  in  the 
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(u2,  v2,  0) 


Figure  5 1 .  Entering  Face  Showing  Node  Coordinates  in  Both  Local  uvw  and 
in  Triangle  Notation. 


triangle  system,  apex  points  were  labeled  according  the  base  length,  b,  the  height,  h,  and 
y.  The  Figure  also  shows  the  corresponding  points  in  the  (uvw)  system.  As  explained 
earlier,  on  the  entering  face,  u=u(v,  w)  so  that  rotation  from  the  (uvw)  to  the  (x',y')  system 
is  given  by 


/  in  in 

X  =ClxV  +  C2xW 

/ 

y  =CjyW. 


(140) 


161 


Then, 


l  in 

b=ClxV2 

i  in 

h  =  clyw 3 


y  = 


in  in 

C1*V3  +  C2xW3 


Clx  V2 


(141) 


and 


Xtn  =  P/]>3  +  Px(Ci>3  +  C£W3) 


iW*-P,<£v 


(142) 


Since  it  was  carried  over  from  the  triangle  code,  the  two-dimensional  root  solve  routine  is 
written  in  terms  of  b,  h,  and  y.  This  resulting  system  of  equations  is  identical  to  Equations 
(138)  and  (139),  so  the  discussion  is  valid  for  either  choice  of  notation.  The  notation  in 
terms  of  the  uvw  system  is  used  to  derive  the  initial  guess  equations. 

B.2.3.  Derivation  of  2D  System  Using  One-Dimensional  Moments  Functions 
Applying  the  definition  for  recursion  on  number  of  arguments  to  Equations  (138) 
and  (139)  gives 


and 


K  (4)-m,(*J 
MoVJ-KW 


=  0, 


(143) 
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(144) 


/2(V,«)  = 


K(~rj-K(xin-rj 


The  functions  fl  and  £2  are  now  given  completely  in  terms  of  single-argument  functions. 
Further,  the  Jacobian  can  also  be  written  in  terms  of  single  argument  functions. 

The  number  of  two-dimensional  root  solves  required  per  iteration  for  a  given 
problem  is  just  the  average  number  of  input  faces  per  tetrahedron  times  the  number  of 
tetrahedron  in  the  problem  times  the  number  of  angles  in  the  quadrature.  In  general,  all 
tetrahedra  are  either  case  ±3  or  4,  so  that  the  average  number  of  input  faces  is  two.  For 
an  S8  quadrature,  there  are  80  angles,  so  that  there  are  160-times-the-number-of- 
tetrahedra  two-dimensional  root  solves  per  iteration.  For  the  system  in  the  form  of 
Equations  (138)  and  (139),  each  Newton  iteration  requires  2  two-argument  moment 
function  evaluations  and  4  three-argument  moment  function  evaluations  to  compute  the 
function  and  jacobian  values.  In  the  single  argument  form,  a  Newton  iteration  requires  4 
single-argument  moment  function  evaluations.  This  is  much  more  computationally 
efficient  than  the  multi-dimensional  formulation.  The  equations  are  numerically  well 
conditioned  provided  the  values  of  and  Yin  are  far  enough  apart  in  magnitude. 
Empirical  testing  shows  that  if  [  Xin  -  4I>  10~4  the  single-argument  formulation  gives 
as  much  precision  as  the  multi-dimensional  equations. 
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B.2.4.  Initial  Guess  Equations  for  the  2D  System 

Although  Newton's  method  is  quadratically  convergent,  if  the  initial  guess  is  not 


sufficiently  close  to  the  root,  the  method  may  take  a  large  number  of  iterations  to 
converge,  and  may  not  even  converge  at  all.  The  initial  guess  generator  ensures  that  the 
two-dimensional  root  solving  routine  converges  in  a  minimum  number  of  iterations  by 
generating  an  initial  guess  for  and  Ym  that  is  close  to  the  correct  value.  To  understand 
how  this  works,  we  must  first  examine  the  ranges  of  py,  p^,  and  their  sum,  psum.  The 
equations  defining  spatial  moments  on  tetrahedral  faces  are 


(146) 


Using  Equation  (146), 
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1 


p,= 


2  Jdwjdv(l-w)  ^;"[v(v?w),^(v7w)] 
0  0 


1  w 

2  Jdw  Jdv  t|F" [v(v,vt>),w(v,vi>)] 

o  o 


0<P ><1 


(148) 


Using  Equations  (137)  and  (140), 


V 


V2 


c 

I  m 

”2 '"3 


(149) 


Again,  using  Equation  (146), 


1  w 


2  I  dw  I  dv(w-v)  i|f'"[v(v,w>),w'(v,v£’)] 

p  =  — £ - 2 - -  0  <  p  <1 

r*y  i  w  Kjcv 

2  (dw  (dv  ijf'” [v(v,m)),w(v,vP)] 


hr 

0  0 


xy 


(150) 


Finally,  define  psum  =  py  +  p^  ,  and 


1  W 

2  Jdw  J dv(l -v)  i|i'”[v(v,w),w(v,w)] 
o  o 


1  W 

2  Jdw  Jdvifr'n[v(v,w),w(v,w)] 

o  o 


0  <  P  sum  <  1  ‘  (151) 
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Equations  (148),  (150),  and  (151)  define  the  area  in  (px,  p^)  space  which  bounds  all  the 
values  of  px  and  pX7  that  could  satisfy  the  2D  system.  This  is  shown  in  Figure  52. 


Figure  52.  Area  in  (py,  pX7)  Space  Containing  Valid  Solutions  to  the  2D 
Rootsolve  System. 

Note  that  the  area  of  possible  solutions  for  Newton's  method  over  a  triangular  face  is  itself 
a  triangle. 

The  general  form  of  the  initial  guess  equations  derives  from  an  examination  of  the 
behavior  of  the  solution  as  the  points  (py,  p^)  move  toward  the  comers  of  the  triangle. 

As  both  py  and  pX7  approach  zero,  the  values  of  and  Ym  that  satisfy  Equations  (138) 
and  (139)  are  given  by 
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-4  ~  ,  and  , 


(152) 


If  both  p  values  are  less  than  0.06,  the  above  equations  give  and  values  that  agree 
with  the  exact  values  to  12  digits,  the  current  tolerance  for  the  root-solved  values.  As  the 
value  of  py  approaches  1, 


X, 


1  “Pi 


,  and  , 


(153) 


If  the  value  of  py  is  greater  than  0.95,  the  above  equations  give  both  X  and  Y  to  12  digits. 
For  the  case  where  p^,  approaches  1, 


4* 


l“Pi 


4 


(154) 


For  values  of  pX7  greater  than  0.95,  the  above  equations  give  X^  and  Ym  to  12  digits  of 
accuracy.  General  case  equations,  that  give  good  approximations  over  the  entire  area, 
can  be  interpolated  from  Equations  (152),  (153),  and  (154)  by  assuming  those  equations 
represent  limiting  cases.  The  general  case  equations  are  then 
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(155) 


yinitial  O  “pL)"  ~  p],? 

= - - - - —  ,  and  , 


1  -P, 


y 


Y: 


initial 


(l-p 2yf 


xy 


When  the  rootsolver  is  called,  it  computes  the  p  values  from  the  input  flux  moments.  If 
the  p's  are  such  that  Equations  (152),  (153),  or  (154)  apply,  those  equations  are  used  to 
generate  and  Ym  without  root-solving.  Otherwise,  Equations  (154)  are  used  to 
generate  an  initial  guess  for  Newton's  method.  Using  initial  guesses  generated  by 
Equations  (154),  the  root-solver  is  generally  able  to  find  and  Y-m  values  good  to  10 
digits  using  fewer  than  four  iterations. 

B.2.5.  Computation  of  Distribution  Coefficients  from  Root-Solved  Values 

Given  the  root-solved  constants,  Px  and  py  are  found  by  inverting  Equations  (133) 


P,= 


x*  -  4 


0-Y4/V4 

h 


(156) 


and 


a  = 


(157) 


168 


Equations  (156)  work  well  regardless  of  the  values  of  and  Ym,  the  same  is  not  true  for 
Equation  (157).  For  large  negative  values  of  and  Ym,  numerical  problems  caused  by 
over  and  underflows  can  produce  results  that  are  not  machine  representable.  For  example, 
if  X^  becomes  a  large  negative  number,  regardless  of  the  value  of  Y^,  the  exponential 
term  could  evaluate  as  zero  or  an  underflow,  both  of  which  cause  the  calculation  to  crash. 
In  this  cases,  an  alternate  form  of  Equation  (157)  is  used.  For  example,  if  X^  =  -750  and 
Ym  =  50,  then  e'V"  =  1.90><10~326,  an  underflow,  M0(Xjn,Yjn)  =  8.76*  103]9,  an  overflow, 
but  the  product  is  1.6666666666666667x1  O'6,  a  result  that  could  be  significant  in  the 
solution  of  the  problem. 

Using  the  definition  of  the  exponential  moments  functions. 


For  large  negative  values  of  X^,  all  but  the  first  term  on  the  right  can  be  neglected,  and 
Equation  (158)  can  be  approximated  as 


'KixinJin) 


XM-XJ 


(159) 


Empirical  testing  shows  that  the  approximation  can  be  used  for  X^  <  -30  and  X^  -  Ym  <  - 
30  with  a  maximum  relative  error  of  10'12  (the  tolerance  used  to  compute  the  moments 
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functions).  Using  the  example  values  of  =  -750  and  Ym  =  50,  Equation  (159)  gives 
the  correct  result  of  1.6666666666666667xl0‘6 ,  with  the  added  benefit  of  not  having  to 
compute  an  exponential  moment  function.  Therefore,  Equation  (159)  is  used  whenever 
possible  to  compute  a. 

If  the  value  of  is  a  large  negative  number  and  the  conditions  for  using  Equation 
(159)  don't  hold.  Equation  (157)  can  be  modified  using  the  identity 

=  •  (iso) 

For  example,  if  X^-750  and  Yjn=-751,  then  e  X'n  =  1.90xl0“326,  an  underflow,  and 
M0(Xjn,Yjn)  =  1.20xl0323  ,  an  overflow.  The  product  of  the  two  is  2.28621637167xl0‘3, 
a  perfectly  reasonable  result.  If  we  instead  use  Equation  (160) ,  ex"  "7”  =  2.72, 
M0(Xw-Ytn, -YJ  =  8.41xl0‘4,  and  the  product  of  the  two  is  2.28621637167xl0'3 , 
thereby  avoiding  overflow/underflow  problems. 
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B.3.  Calculation  of  the  Source  Distribution  Coefficients:  3D  System 

B.3.1.  Derivation  of  3D  System  Using  Multi-Dimensional  Moment  Functions 

The  assumed  form  of  the  source  distribution  for  the  EC  method  is 


S(x,y,z)  =  aebx'cy  +  dz 


(161) 


For  a  given  tetrahedron,  the  source  moments  are  calculated  using  the  scalar  fluxes  from 
the  previous  iteration.  The  source  distribution  coefficients  are  then  chosen  to  match  the 
zeroth  and  first  source  moments.  Using  the  knowledge  of  the  source  distribution 
moments,  the  coefficients  of  the  distribution  can  be  computed  by  inverting  Equations  (84), 
giving  a  system  of  three  equations  in  three  unknowns  which  then  has  to  be  rootsolved. 

The  equations  are 


Pw  = 


K  (XS,YS,ZS) 


M,,SA  M0(Xs,Ys,Zs) 


(162) 


v^A  v2  w  M0(YS-XS,ZS-XS,-XS) 


(163) 


Pu  = 


Su 

u1S' 


U3V2~U2V3 


uxv2 


PW  +  — Pv 


U , 


ms-zs,xs-zs,-z) 


(164) 


171 


Now  define 


P  w 


Mx(yrx„zrx,t-x) 

M0(YS-XS,ZS-XS,-XS) 


(165) 


K(Y~Zs>X-Zs>-Zs) 


(166) 


giving  the  system  of  equations 


fr(Vs,Zs)  =  pM 


K  &,>y„z,)_o 

M0(XS,YS,ZS ) 


(167) 


fJX,Ys,Zs)  =  pv 


m0(ys-xs,zs-xs,-xs) 


(168) 


and 


fXXYZs)  =  o 


^Xs^r^-z) 

K(Ys~zs,xs-zs,-zs: ) 


(169) 


The  above  system  contains  three  equations  in  the  three  unknowns  Xs,  Ys,  and  Zs. 
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B.3.2.  Derivation  of  3D  System  Using  One-Dimensional  Moment  Functions 

Just  as  in  the  two-dimensional  case,  the  three-dimensional  system  can  be  written 
entirely  in  terms  of  single-argument  function  using  recursion  on  number  of  arguments. 
This  leads  to  the  following  system  of  equations: 


f\ =  Pw  ~ 


(Y-Z)M0(Xs)  +  (Zs-Xs)M0(Ys)+(Xs-Ys)M0(.Zs) 


=  0,  (170) 


fJX,Y,Zs)  =  Pv 


(Z,  -  W  (-*,)  -  W,  -X,) +  Y_M&  -*,) 

(Z,  -  YSW0(-XS)  -  ZSM0(YS  -xs)  H-  YSM0(ZS  -Xs) 


(171) 


and 


UXsJZ)  =  9u 


YMi&rZj  -xW'-z,)  +  (r,-r,  w^j-zj 

YM0(Xs-Zs)-XsM0(Ys-Zs)  +  (Xs~Ys)M0(-Zs) 


(172) 


Now  both  the  functions  and  their  derivatives  used  in  the  Jacobian  are  given  completely  in 
terms  of  single-argument  moment  functions.  One  Newton  iteration  using  the  multi¬ 
dimensional  form  of  the  system  requires  three  calls  to  the  three-dimensional  routines  and 
nine  calls  to  the  four-dimensional  routines.  In  the  single-argument  form,  a  Newton 
iteration  requires  only  nine  calls  to  the  single-argument  routines.  As  in  the  2D  case,  the 
above  system  is  numerically  ill  conditioned  as  any  of  the  variables  Xs,  Ys,  Zs,  approach 
zero,  or  as  the  difference  between  any  two  variables  approaches  zero.  The  system  in 
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terms  of  multiple-argument  functions  agrees  with  the  multi-dimensional  formulation  when 
any  of  the  arguments  has  magnitude  less  than  1CT4,  or  when  the  difference  between  any 
two  arguments  is  less  than  KT4  .  Since  it  is  much  more  efficient,  the  system  in  terms  of 
single-argument  functions  is  used  whenever  possible. 

B.3.3.  Initial  Guess  Equations  for  the  3D  System 

The  similarity  between  the  two-  and  three-dimensional  systems  suggests  that  the 
initial  guess  equations  would  share  some  of  the  same  properties.  This  is  indeed  the  case. 
The  equations  for  generating  an  initial  guess  are  found  using  the  same  procedure  as  the  2D 
case,  beginning  with  an  examination  of  the  ranges  of  pw,  pTO,  and  puvw.  The  equations 
defining  the  source  moments  over  the  tetrahedron  volume  are 

v2<»  “oJv.w) 

dw  f  dv  f  du 

J  v,  J  u, 

Vj(w)  2  Um(v,w)  1 

Using  several  changes  of  variables. 


w3  v2<» 

.{dw  f  dv  f 

J  W,  J  V,  J 

0  i  v,(w)  2  !<„(v,w) 


S(u,v,w) 


1  W  V 


(174) 


6  jdw  J'dvj'du  S[u(u,v,w),v(u,v,w),w(u,v,w)] , 


0  0  0 


where 
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w(u,v,w)  =  (1  -m>)w3 

v(u,v,w)  =  v3  -v2v  +  (v2+v3)m>  (175) 

u(u,V,w)  =  M3  +MjM  -  M2V  +  (tt2  -w3)w  . 


From  Equation  (162), 


Pw  = 


^  c1  A 

w3S A 


1  W  V 

6'j'dw  j'dvj'du  (1  -w)  S(u,v,w) 

0  0  0 _ 

1  W  V 

6jdwjdvjdu  S(u,v,w ) 

0  0  0 


o<P  <1 


(176) 


From  Equation  (165), 


1  W  V 

6  fdw  f  dv  fdu  (v-w)  S(«,v,w) 


„ _  3  _  _  0  0  0 

Pvw  Pv  Pw  1 

v  1  W  V 

V2 


0<P™<1 


6j  dw  j  dv Jdu  S(u,v,w) 
0  0  0 


From  Equation  (166), 


(177) 


1  W  V 


u2 

Pu~  —  Pv- 

UX 


•••  0<P„vW<1(178) 


0  0  0 
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Finally,  the  sum  of  the  p's. 


p  =p+p  +  p  = 

*  sum  •  w  *  vw  "  uvw 


1  W  V 

ij'dwj'  dvj i 
0  0  0 


du  (1  -v  +  u)S(u,v,w) 


1  W  V 

i  J  dw  j  dv  j't 

0  0  0 


•••0<P^<1-  (179) 


du  S(u,v,w ) 


These  equations  define  the  volume  in  pw,  pw  puvw  space,  shown  in  Figure  53,  that 
contains  all  values  (pw,  pw  puvw)  that  satisfy  the  system.  Note  that  in  analogy  with  the 
two-dimensional  case,  the  volume  of  attainable  values  for  the  p's  of  the  source  distribution 
in  a  tetrahedral  cell  is  itself  a  tetrahedron. 

The  general  equation  for  initial  values  is  found  by  examining  the  behavior  of  the 
system  as  the  points  (pw,  pTO!  puvw)  move  to  the  tetrahedron  apexes.  As  the  p's  approach 
zero,  the  following  equations  hold 


and 


(181) 


For  p^  pw  and  puvw  all  less  than  4.0xl0'2,  the  Equations  (180)  give  Xs,  Ys,  and  Zs  to 
twelve  digits,  and  rootsolving  is  not  necessary.  As  pw  approaches  1, 
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(181) 


and  , 


P  uvw 


For  Pw  >  0-95  and  psum  >  0.99,  Equations  (181)  gives  Xs,  Ys,  and  Zs  to  twelve  digits,  and 
rootsolving  is  not  necessary.  As  pTO  approaches  1, 


and  , 


1 


1-p 


sum 


l  P 

UVW  I  w 


(182) 


For  Pvw  >  9  95  and  psum  >  0.99,  Equations  (182)  gives  Xs,  Ys,  and  Zs  to  twelve  digits,  and 
rootsolving  is  not  necessary.  As  puvw  approaches  1, 


and  , 


1 


1-p 


sum 


(183) 


For  puvw  >  0.95  and  psum  >  0.99,  Equations  (183)  gives  Xs,  Ys,  and  Zs  to  twelve  digits, 
and  rootsolving  is  not  necessary.  Assuming  that  the  behaviors  in  the  comers  are  limiting 
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The  initial  values  generated  by  Equations  (184)  generally  allow  the  rootsolver  to  converge 
in  fewer  than  four  iterations. 


Figure  53.  Phase  Space  Volume  Containing  Solutions  to  3D 
RootSolving  System. 


178 


B.3.4.  Computation  of  Distribution  Coefficients  from  Rootsolved  Values 

Given  the  root-solved  values  of  Xs,  Ys,  and  Zs,  the  values  of  the  exponential 
distribution  coefficients  are  found  by  inverting  Equations  (84),  giving. 


Y-Z 

b  =  — - - 

u , 


-uxXs  +  («,  -u2)  Ys  +  u2Zs 


Uiv2 


(185) 


d  = 


_  Ul  V3  Xs  +  (Ul  V2  '  M3  V2  ~  Ul  V3  +U2Vz)Ys  +  (W3  V2  “  U2  V3  )Zs 


uxv2w3 


and , 


a  = 


6er’M0(Xs,Ys,Zs) 


(186) 


For  extreme  values  of  Xs,  Ys,  and  Zs,  Equation  (186)  is  inadequate.  If  Ys  is  a  large 
negative  number,  the  exponential  term  could  be  a  very  small  number,  while  the  moments 
function  is  a  very  large  number.  Although  when  multiplied  together  the  product  of  the 
two  is  reasonable,  the  machine  may  assign  the  exponential  term  either  a  value  of  zero,  or 
underflow,  either  of  which  would  produce  a  non-machine-representable  result.  To  avoid 
this,  alternate  forms  of  Equation  (186)  are  used  based  on  the  values  of  Xs,  Ys,  and  Zs. 
Applying  the  definition  of  the  exponential  moment  functions  and  expanding. 
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(187) 


er’M0(X,Ys,Zs)  = 


Ys  Y-Xs 

e  -  e 


1  -  e ' 


X, ( r, -X,)(Z,  -X)  YS(YS - X ,)  (Zs - X ,) 


1  -  £ J 


r,  y '-zs 

e~e 


Y,V,-Y,)V,-V  Z,(Z,-YJ(Z,-XS) 


If  Ys  <  -30,  Ys-Xs  <  -32,  and  Ys-Zs  <  -32,  then  Equation  (187)  can  be  approximated  to 
good  precision  using 


e**M0{XsJs,Zs)» 


1 


1 


Y,{YrXt){ZrX)  Ys (Zs -  Ys)  (Zs -Xs) 

1 


Y,(YM-XtHZrY) 


(188) 


For  example,  if  Ys=  -750,  Xs  =  -50  and  Zs  =  -51,  then  eYs  =  1.901*lO~309,  which  could 
be  interpreted  as  either  zero,  or  an  underflow,  either  of  which  blows  up  Equation  (186). 
However,  for  these  values,  e¥sM0(Xs,Ys,Zs)  =  2.72498 126*10~9.  Using  the 
approximation.  Equation  (188),  gives  the  correct  result,  2.72498126xl0'9’  but  avoids 
non-numerical  results.  Further,  Equation  (188)  doesn't  require  the  computation  of  a  three- 
argument  moment  function,  which  is  much  faster.  Therefore,  Equation  (188)  is  used 
instead  of  Equation  (186)  whenever  possible. 

For  situations  where  all  the  arguments  of  the  moments  function  in  Equation  (186) 
are  much  less  than  zero,  but  the  conditions  for  using  Equation  (188)  don't  hold,  the 
identity 


YS,ZS) 


eY*  X‘M0(YS 


-XZ  -X  ,-X) 

S 7  S  S9  S' 


(189) 
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can  be  used.  For  example,  if  Xs  =  -1001,  Ys  =  -1000,  and  Zs  =  -1002,  e¥*=  5. 075x1  O'435, 
and  M0(Xs,Ys,Zs )  =  2.899xl0431,  but  the  product  of  the  two  is  0.00147158747978. 
Equation  (189)  becomes  e  1M0(  1,2, 1001)  =  0.00147158747978,  which  gives  the  correct 
result  without  producing  overflows  or  underflows. 


Appendix  C.  Tetrahedral  Splitting 


The  splitting  of  arbitrarily  oriented  tetrahedra  along  the  streaming  direction  can  be 
broken  down  into  two  steps.  The  first  step  is  case-dependent  and  consists  of  determining 
the  identification  of  the  nodes  (vertices)  defining  each  sub-tetrahedron.  The  second  step  is 
case-independent  and  consists  of  constructing  the  rotation  and  translation  matrices  needed 
to  rotate/translate  the  node  coordinates  to  the  local  sub-tetrahedron  coordinate  system, 
computing  the  subcell  coordinates  in  the  subcell  local  coordinate  system,  and  computing 
the  sub-tetrahedron  entering  and  exiting  face  areas  and  volume.  Since  the  splitting 
routines  are  designed  to  support  the  second  step,  this  section  begins  by  discussing  the 
case-independent  second  step  first,  and  then  examines  the  details  of  determining  the 
rotation  matrices  and  sub-tetrahedron  node  coordinates  for  each  of  the  four  cases. 

C.l  Construction  of  the  Local  Sub-Tetrahedron  Coordinate  System 

The  splitting  routines  identify  sub-tetrahedra  that  are  arbitrarily  oriented  with 
respect  to  the  centroid  coordinate  system,  an  example  of  which  is  shown  on  the  left-hand 
side  of  Figure  54.  Application  of  the  quadrature  requires  the  construction  of  a  local  sub¬ 
tetrahedron  coordinate  system  as  shown  in  the  right-hand  side  of  the  figure.  For  each  sub- 
tetrahedron,  the  splitting  routine  generates  a  node  matrix,  T^,  that  contains  the  centroid- 
system  coordinates  of  nodes  defining  the  tetrahedron. 


T  = 


' A  B  C 

x  X  X  X 

Ay  By  Cy  D, 
B,  Cz  DzJ 


(190) 
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Figure  54.  Case  1  Tetrahedron  in  Global  xyz  Coordinate  System  (left  figure)  and  in  Local 
uvw  Coordinate  System  (right  figure) 


The  nodes  are  specifically  ordered  so  that  when  transformed,  the  local  coordinate  matrix 
has  the  form 


T  = 

UVW 


{n  \ 
0  Wj  u2  u3 

0  0  v2  v3 


0  0  0  w. 


V 


(191) 


The  transformation  is  accomplished  using 


(192) 
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where  R,  the  rotation  matrix  from  the  centroid  coordinate  system  to  the  local  uvw 
coordinate  system,  is  defined  as 


R  = 


e  me  e  ~e  e  -e 

U  X  V  X  W  X 


e  -e  e  -e  e  -e 

u  y  v  y  w  y 


f  £  £  *  /\  /\  /V 

[  e  -e  e  *e  e  -e 

\  u  z  V  z  w  z 


and  the  translation  vector 


t  =  R 


'A.' 


(193) 


(194) 


is  just  the  origin  of  the  local  coordinate  system  in  rotated  coordinates.  (Note  that  the 
translation  is  accomplished  after  the  rotation  because  for  case  4  tetrahedra  the  translation 
vector  is  unknown  until  after  the  rotation).  The  rotation  matrix  is  defined  in  terms  of  dot 
products  to  avoid  repeated  computation  of  sines  and  cosines. 

The  rotation  matrix  is  composed  of  dot  products  between  the  unit  vectors  along 
the  coordinate  axes  of  each  system.  Since  the  centroid  coordinate  system  is  simply  the 
global  coordinate  system  translated  to  the  centroid,  the  unit  vectors  are  the  same  as  those 
of  the  global  system.  In  the  local  uvw  system,  the  u-axis  is  defined  by  the  streaming 
direction.  The  unit  vector  in  the  u-direction  is  simply 

eu=Q  =  \iex  +  r] ey  +  £ e,  .  (195) 
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Referring  to  Figure  54,  since  the  point  C  lies  in  the  uv-plane,  the  vector  PAC  must  also  lie 
in  the  uv-plane  and  the  unit  vector  in  the  w-direction  can  be  defined  as 


e  =e  x 

W  U 


AC 


\P 


AC  l 


(196) 


then 


ev  =  e*Xeu-  (197) 

R  and  t  are  also  used  to  transform  the  source  distribution  coefficients  to  the  local  system, 
and  to  transform  the  cell  flux  moments  back  to  the  centroid  system. 

C.2  Case  4  Tetrahedra  Splitting 

Case  4  tetrahedra  are  characterized  by  having  two  input  and  two  output  faces. 
When  split  along  the  streaming  direction,  the  case  4  tetrahedron  produces  four  case  1 
tetrahedra  as  shown  in  Figure  55. 

In  essence,  the  problem  of  splitting  a  Case  4  tetrahedron  reduces  to  finding  the 
coordinates  of  the  points  P5  and  P6  in  Figure  56.  Given  these  points,  the  matrices 
containing  the  nodes  and  the  rotation  and  translation  matrices  can  be  generated.  There 
are  several  ways  to  approach  the  problem  of  finding  points  P5  and  P6.  One  approach  is  to 
note  the  following: 

1)  Point  P5  is  on  the  line  P,P3 

2)  Point  P6  is  on  the  line  PJt>4 
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3)  Point  P5  is  on  the  line  parallel  to  Q  passing  through  P6 
and  solve  the  resulting  system  of  six  equations  for  the  coordinates  of  P5  and  P6.  Another 
approach  is  to  note: 

1)  Point  P5  lies  in  the  plane  defined  by  P^Ps  and  the  plane  defined  by  P!P3P4 

2)  Point  P6  lies  in  the  plane  defined  by  P1P2P4  and  the  plane  defined  by  P2P3P4 

3)  Point  P5  is  on  the  line  parallel  to  Q  passing  through  Point  P6. 

and  solve  the  resulting  set  of  equations  for  the  P5  and  P6  coordinates.  Both  of  these 
approaches  were  investigated,  in  each  case  the  resulting  equations  for  the  node 
coordinates  were  very  large  and  unwieldy,  and  did  not  always  work. 


Figure  55.  Case  4  Tetrahedron  and  Resulting  Sub-Tetrahedra. 
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Figure  56.  Case  4  Tetrahedron  Showing  Point  Labeling  Convention 
(Configuration  1). 


The  approach  that  finally  did  work  relies  on  the  fact  that  since  the  translation. 
Equation  (194),  is  defined  in  terms  of  rotated  coordinates,  then  it  is  only  necessary  to 
know  P5  and  P6  in  the  rotated  coordinates,  and  the  rotation  can  be  done  first.  Once 
rotated,  the  problem  of  finding  the  intersection  of  three  lines,  P{P3 ,  P2P4  and  a  line 
parallel  to  the  streaming  direction,  in  three  dimensions  can  be  reduced  to  finding  the 
intersection  of  two  lines  in  a  plane. 

The  first  step  is  to  label  the  tetrahedron  nodes.  When  passed  to  the  splitting 
routine,  the  coordinates  of  the  parent  tetrahedron  nodes  are  known,  but  we  have  no  idea 
how  the  nodes  are  oriented  with  respect  to  one  another.  The  nodes  are  labeled  using  the 
following  convention,  illustrated  in  Figure  56.  The  point  where  the  streaming  vector 
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enters  the  tetrahedron  is  labeled  P5,  and  the  exiting  point,  P6.  The  nodes  defining  the  edge 
containing  P5  (the  edge  common  to  the  two  entering  faces)  are  labeled  Pj  and  P3.  The 
nodes  defining  the  edge  containing  P6  (the  edge  common  to  the  two  exiting  faces)  are 
labeled  P2  and  P4.  There  are  two  possible  point  configurations  which  can  result:  the  one 
shown  in  Figure  56,  and  one  similar  to  Figure  56  in  which  the  points  Pj  and  P3  are 
exchanged.  (Another  approach  is  simply  to  swap  labels  Pj  and  P3,  the  method  described 
is  slightly  more  efficient  and  straightforward).  The  routine  first  checks  the  orientation  and 
uses  a  set  of  equations  consistent  with  the  orientation.  For  configuration  1,  shown  in 


Figure  56,  the  sub-tetrahedra  are:  T^PjTWPo),  T2=(P5,P6,P2,P3),  T3=(P5,P6,P3,P4), 


and  T4-(P5,P6,P4,P1).  And  for  configuration  2  the  sub-tetrahedra  are:  T1=(P5,P6,P2,P1), 
T2=(p5>p6>p3>p2)>  T3=(p5,p6,p4,p3), T4=(P5,p6,p1  ,P4 ).  Again,  we  note  that  the  node 
ordering  is  important  in  identifying  the  positions  of  the  nodes  once  transformed  into  the 
subcell  system. 

C.2.1  Sub-Tetrahedron  1  (Configuration  1) 

Before  transforming  coordinates,  we  recognize  that  sub-tetrahedron  1  is  defined  as 

(*S  *6  *1  ^ 


>5  t6  yl  y2 


(198) 


‘2) 


where  the  node  coordinates  are  in  the  centroid  coordinate  system.  The  node  coordinates 
are  given  in  the  centroid  coordinate  system  because  the  resulting  rotation  and  translation 
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matrices  will  also  be  used  to  transform  the  source  moments  from  the  centroid  system  and 
to  transform  the  subcell  moments  back  into  the  centroid  system.  A  typical  example  is 
shown  in  Figure  57,  where  the  xyz  axes  represent  the  centroid  coordinate  system,  the 
u'v'w'  axes  represent  the  rotated  (but  not  yet  translated)  coordinate  system,  and  the  uvw 
axes  the  folly  transformed,  local  sub-tetrahedron  coordinate  system. 


Figure  57.  Case  4,  Sub-Tetrahedron  1  in  Relation  to  Global,  Rotated, 
and  Local  Sub-Tetrahedron  Coordinate  Systems. 


For  sub-tetrahedron  1  (in  fact  for  all  four  sub-tetrahedra),  by  construction: 

eu=Q  =  iiex  +  r\ey  +  ^ez.  (199) 


Further,  referring  to  Figure  58,  Pj  and  P3  both  lie  in  the  plane  perpendicular  to  ew,  thus 
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(200) 


<?...  =  Qx 


P, 


31 


I P 


31  I 


whereP31  is  the  vector  from  point  P3  to  point  Pj  and 


e  =  e  x  e  . 

V  W  U 


(201) 


(Note  that  eu  =  eui,  ev  =  ev,,  and  ew  =  ew,.)  We  can  generate  the  rotation  matrix  using  the 
unit  vectors  in  Equation  (193)  and  rotate  the  entire  parent  tetrahedron  into  the  u'v'w' 
system,  as  shown  in  Figure  58.  The  node  coordinates  are  now  given  in  the  rotated 
coordinate  system. 

Next,  examine  the  projection  of  the  edges  into  the  v'w'-plane,  shown  in  Figure  59. 
From  the  line  connecting  P'j  and  P'3,  it  is  obvious  that, 

W5  =w6  =w\  =w/3  -  (202) 


The  equation  of  the  line  from  P'2  to  P'4  is  just 


w 


/ 


(w2  -w'4) 
(v'3-v'4) 


(v' -v\)  +w'4  . 


(203) 


Using  Equation  (202)  and  applying  Equation  (203)  at  the  point  P’5  (i.e.  let  w  ’  =  w  '5 ,  then 
v '  =  v  75 ,  and  solve  for  v  '5 ) 
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(204) 


v 


/ 

5 


(V i  -M^) 

("VO 


(v'2-v\)  + 


/ 

4  ‘ 


From  the  projection  of  the  edge  connecting  P’j  and  P'3  into  the  u'w'-plane. 


u 


O  3  ~m/i) 

(V3-v/l) 


K-v'). 


(205) 


Using  Equations  (202)  and  (205)  gives. 


u 


/ 

5 


O3  -«;i) 

K-v') 


K-vO 


(206) 


Again,  in  the  u'w'-plane. 


u 


(w1  -w'2) 
O4  -w72) 


(«  4 -“a). 


(207) 


from  which 


/ 

6 


(w'4-w'2) 


(208) 
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Figure  58.  Case  4  Tetrahedron  in  Sub-Tetrahedron  1  Rotated  Coordinate  System. 


Figure  59.  Projection  of  Case  4,  Sub-Tetrahedron  1  Edges  into  the  v'w'-Plane. 
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We  now  have  all  the  elements  of  the  sub-tetrahedron  coordinate  matrix 


1  /  /  /  / 

U  5  »  6  W  1  U  2 


V75  v76  V7x  V72 


Vh,75  <  w,/i  wy 


(209) 


and 


(210) 


At  this  point  all  the  information  is  available  to  transform  the  tetrahedron  into  the  local  uvw 
system. 

Since  the  position  of  the  u-axis  is  the  same  for  all  four  sub-tetrahedra,  the  values  of 
u'5  and  u'6  are  the  same  for  all  sub-tets.  Further,  since  v'5=v'6  and  w'5=w'6,  we  only  need 
compute  v'5  and  w'5  and  the  unit  vectors  for  each  of  the  remaining  sub-tets  to  complete 
the  splitting  procedure. 

C.2.2 .  Sub-Tetrahedron  2  (Configuration  1) 

For  sub-tetrahedron  2, 


~  ~  x  42 

e  ~  e  x  — — 

w  u  .  .  5 


and 


1^ 


e  -  e  x  e 

V  w  u 


(211) 
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Again,  rotate  the  entire  parent  tetrahedron  into  the  new  u'v'w'  coordinate  system.  From 


the  projection  of  P{P3  and  P2P4  into  the  v'w'  plane,  w/2  =  w'4  =  w  '5=<  ,  and 


v<  = 


(w'  -  W73) 

K-<) 


(v/3  -V\)+V‘ 


(212) 


which,  at  P’5,  gives 


v 


/ 

5 


(V5  -^3) 


~V\)  + 


3  • 


(213) 


C.2.3.  Sub-Tetrahedron  3  (Configuration  1) 

For  sub-tetrahedron  3, 


5  A  A  A  „  13 

e  =  U ,  e  =  e  x  — — 

«  ’  w  u  .  ^  ’ 


and 


I* 


e  -  e  *  e 

V  w  w 


13  I 


(214) 


By  symmetry,  w'5  and  v'5  can  be  found  using  w\=w'^=wl5=w!6  and  Equation  (204). 

C.2.4.  Sub-Tetrahedron  4  (Configuration  1) 

For  sub-tetrahedron  4, 


*  ~  24 

e  —  e  x 
w  a  .  -  .  > 


and 


I* 


e  -  e  x  e 

V  w  u 


(215) 
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By  symmetry,  w'5  and  v'5  can  be  found  using  w'2=w'4  =  wl5  =  w,6  and  Equation  (213). 
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C.2.5  Configuration  2 

The  equations  for  Configuration  2  can  be  derived  in  a  completely  analogous 
manner  to  the  Configuration  1  equations.  The  resulting  equations  are  identical  to  the 
Configuration  1  equations  with  the  exception  that  u'5  and  u'6  are  computed  on  Sub- 
Tetrahedron  1  using 


u'5=u'  j  + 


(w75  -w;i) 


(216) 


and , 


»6=«2 


(V5~V2) 

K-vg 


(u'-u'2) 


(217) 


C.3.  Case  3  (-3)  Tetrahedra  Splitting 

Case  3  and  -3  tetrahedra  are  characterized  by  having  either  three  input  and  one 
output,  or  one  input  and  three  output,  faces.  When  split  along  the  streaming  direction,  the 
Case  3  (-3)  tetrahedron  subdivides  into  three  Case  1  sub-tetrahedra,  as  shown  in  (217).  A 
Case  2  (-2)  tetrahedra  results  when  the  streaming  direction  intersects  the  entering  (exiting) 
face  along  one  of  the  edges,  and  a  Case  1  tetrahedron  results  when  the  intersecting  point  is 
coincident  with  one  of  the  entering  (exiting)  face  apexes. 
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Figure  60.  Case  3  Tetrahedron  and  Resulting  Sub-Tetrahedra 


Each  tetrahedron  is  defined  by  four  nodes,  each  having  an  x,  y,  and  z  coordinate  in 
the  global  xyz  coordinate  system.  There  is  no  a-priori  way  of  knowing  how  the  individual 
nodes  are  oriented  with  respect  to  one  another  in  space,  so  the  following  labeling 
convention  is  used.  When  a  tetrahedron  is  passed  to  the  Case  3  (-3)  routines,  the  nodes 
defining  the  entering  (exiting)  face  are  labeled  P1?  P2,  and  P3,  and  the  remaining  node  is 
labeled  P4  P4  is  identified  by  simply  summing  the  four  node  numbers  defining  the 
tetrahedron,  and  subtracting  the  three  node  numbers  defining  the  entering  face.  The  result 
is  the  node  number  of  P4.  Once  the  intersection  point  between  the  streaming  direction  and 
the  entering  (exiting)  face  is  determined,  it  is  labeled  P5.  This  results  in  four  possible 
configurations,  shown  in  Figure  61,  where  the  configuration  numbers  refer  to  the 
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orientation  of  points  P]3  P2,  and  P-  about  the  point  P5  on  the  outward  side  of  the  entering 
(exiting)  face.  The  orientation  is  determined  using 

\  PU  x  Pu  j  •  n  >0,  Configuration  1 
\  Pn  xPl3 )  •  n  <0,  Configuration  2 
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where  P 12  and  P 13  are  vectors  from  Pj  to  P2  and  from  P}  to  P3  respectively,  and  n  is  the 
entering  (exiting)  face  outward  normal.  Once  the  specific  configuration  is  known,  the 
orientation  of  the  tetrahedra  and  the  relative  position  of  all  the  nodes  are  known,  and  the 
sub-tetrahedra  and  rotation/translation  matrices  can  be  generated.  In  the  following 
discussion,  we  use  Case  -3,  Configuration  1  as  an  example.  Figure  62  shows  a  Case  -3 
tetrahedron.  The  first  step  in  splitting  the  tetrahedron  along  the  streaming  direction  is  to 
determine  the  coordinates  of  the  point  P5,  where  the  streaming  direction  intersects  the 
entering  face.  Mathematically  this  is  equivalent  to  finding  the  intersection  between  a  line 
parallel  to  Q  passing  through  point  P4,  and  the  plane  defined  by  points  Pl3  P2,  and  P3.  The 
equation  for  line  parallel  to  Q  passing  through  P4  is 
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(219) 


*5  ~X4  =  T5  ~  y A  =  *5  ~Z4 

P  r)  l 


The  gives  two  independent  equations  in  the  three  unknowns  x5,  y5,  and  z5.  We  also  know 
that  P5  must  fall  in  the  plane  P,  P2  P3,  which  means  that  the  vector  from  P,  to  P5,  Pl5 , 
has  to  be  orthogonal  to  the  face  normal,  n,  where 

n  =  nxex+nyey  +  nzez,  (220) 


and  the  third  equation  is  then 

(x5  -  x,  )nx  +  (y5  -y}  )ny  +  (z5  -  zx )nz  =  0  . 

This  gives  the  following  system  of  three  equations  in  the  three  unknowns  x5,  y5,  and  z5, 

r\(x5  - x4)  -  p(_y5  -y4)  =  0 

Uy5-yA)-r -](z5-z4)  =  o  (222) 

(x5  -  Xj  )nx  +  (y5  -yl)ny  +  (z5  -  zx )nz  =  0  . 

Solving  the  system  in  Mathematica  and  simplifying  using  vector  notation  gives 
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n-P 


*5  P 

y5  =  r\ 

zs  =  Z 


41 


n  •  Q 


+  xA 


n-P 


41 


n  •  Q 

n-P 

n  ■  Q 


>4 


+  Z 


4  • 


(223) 


In  Equations  (223),  g,  q,  x4,  y4,  z4,  and  n  -  Q ,  are  all  known.  The  only  quantity  that 
has  to  be  computed  is  n  •  P4l ,  which  need  only  be  computed  once  to  get  the  three 
coordinates.  This  same  system  of  equations  is  used  for  both  Case  3  and  Case  -3,  and  for 
both  configurations.  Given  the  coordinates  of  P5,  all  the  necessary  information  is  available 
to  generate  the  transformation  matrices  and  transform  the  sub-tetrahedra  into  the  local 
coordinate  system. 

Equations  (223)  can  also  be  obtained  by  geometric  arguments.  Equations  (223) 
are  equivalent  to 


=  Q 


t 


\ 


a  *  p 

n  r  41 

ri-Q  j 


(224) 


Examining  Figure  63  we  see  that  the  altitude  is  b  =  P41  •  n  and  that  the  edge  is 
c-b /  cos0  where  cos0  =  n  -  Q,  thus 


P4l-n 

c  = - -  •  (225) 

n  •  Q 
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Equation  (224)  follows  from  the  fact  that  P45  =  cCl . 


C.4.  Case  2  (-2)  Tetrahedra  Splitting 

Case  2  (-2)  tetrahedra  are  characterized  by  having  either  two  input  and  one  output, 
or  one  input  and  two  output  faces.  When  split  along  the  streaming  direction,  a  Case  2  (-2) 
tetrahedron  produces  two.  Case  1  tetrahedra,  as  shown  in  Figure  64.  The  Case  2 
tetrahedron  is  a  degenerate  case  of  the  Case  3  tetrahedron.  The  labeling  convention  for 
the  Case  2  tetrahedron  is  the  same  as  that  used  for  Case  3,  with  the  exception  that  the 
point  P5  is  restricted  to  lie  along  PXP2 .  As  with  Case  3,  this  results  in  two  possible 
configurations  each  for  Case  2  and  Case  -2.  The  splitting  is  handled  exactly  as  with  the 
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Case  3  tetrahedron,  the  coordinates  of  point  P5  are  found  using  Equation  (223).  Given 
these  coordinates,  all  the  necessary  information  is  available  to  generate  the  transformation 
matrices  and  transform  the  sub-tetrahedra  into  the  local  coordinate  system. 

C.5.  Case  1  Sub-tetrahedra 

Case  1  tetrahedra  are  characterized  by  having  only  one  input  and  one  output  face. 
No  splitting  is  necessary  since  the  tetrahedron  already  has  the  desired  special  case 
characteristics. 
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Appendix  D.  Conservation  Equations 

Solved  for  the  cell  flux  moments,  the  balance  equations  are 


ijr4  = 

i|iw  = 

r  = 

= 


U1  SA+3ip" 

0  ,  out 

_ 

6 
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_  34 C 
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,  in  ~  ,  out 

€ 

e 

(226) 


In  this  form,  the  balance  equations  provide  a  much  more  computationally  efficient  means 
of  computing  the  subcell  flux  moments  since  all  quantities  on  the  right  hand  side  of 
Equations  (226)  are  known.  For  example,  to  compute  the  subcell  flux  moments  for  the 
EC  method  using  Equations  (89)  and  (91)  the  code  must  compute  three  distinct  three- 
argument  moment  functions  and  four,  four-argument  moment  functions.  The  same 
calculation  using  the  balance  equations  is  essentially  free.  Note  that  for  LC,  CLC,  and  SC, 
all  of  the  moments  functions  needed  to  compute  the  subcell  flux  moments  must  be 
computed  to  get  the  source  and  entering  flux  moments,  so  there  is  no  efficiency  gain  in 
using  the  balance  equations  for  these  methods. 

There  is  a  drawback  to  using  Equations  (226).  As  the  cell  optical  thickness 
approaches  zero  (e-0)  they  become  ill-conditioned,  and  loss  of  significant  digits  can 
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occur.  Note  that  the  equations  are  all  in  the  form  a-b.  We  can  estimate  the  precision  lost 
by  examining  the  quantity 


n  =  log 


/ 

i 

a-b  |  ) 

< 

a+b 

i 

(227) 


If  n  is  less  than  zero,  then  loss  of  digits  due  to  subtraction  of  nearly  equal  quantities  has 
occurred  resulting  in  a  loss  of  approximately  n  digits.  We  can  use  Equation  (227)  to  limit 
the  precision  lost  in  computing  the  cell  moments  using  the  balance  equations  by  requiring 

i|ry<  (a  +  b)\0m  (228) 

where  j= A,  u,  v,  or  w,  and  m  is  the  maximum  number  of  digits  allowed  to  be  lost.  For  all 
problems  computed  in  Chapter  5,  the  moments  functions  were  computed  to  twelve  digits, 
and  the  source  and  entering  flux  distribution  coefficients  were  computed  to  ten  digits, 
therefore  it  is  not  reasonable  to  demand  16  digits  of  accuracy  in  the  cell  flux  moments 
computation.  At  most,  the  cell  moments  need  only  be  as  accurate  as  the  moment  function 
computation  tolerance.  Further,  the  final  results  are  less  sensitive  to  the  number  of 
accurate  digits  in  the  cell  first  moments  than  in  the  number  of  accurate  digits  in  the 
average  moments,  hence  the  first  moments  do  not  have  to  be  computed  to  as  much 
accuracy  as  the  zeroth  moments. 

The  EC  quadrature  routine  using  the  balance  equations  works  as  follows: 

Compute  the  exiting  flux  moments 

Compute  the  source  moments  using  the  distribution  coefficients 
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Compute  the  entering  flux  moments  using  the  distribution  coefficients 
Compute  the  cell  flux  moments  using  the  balance  equations 
If  (lost  >  4  digits  for  zeroth  or  >  6  digits  for  first  moments)  then 
Recompute  cell  moments  using  Equations  (89)  and  (91) 

Endif 
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Appendix  E.  Mesh  Generation 


E.l  IDEAS 

Tetrahedral  meshes  were  generated  using  the  IDEAS  (Lawry,  1991)  code.  IDEAS  is  a 
finite  element  code  with  the  capability  to  generate  three-dimensional  meshes  using  a 
variety  of  cell  geometries.  The  IDEAS  Student  Guide  (Lawry,  1991)  and  the  set  of 
IDEAS  manuals  contain  information  on  generating  meshes.  Once  generated,  the  data  are 
written  to  a  file  in  NASTRAN  format.  The  NASTRAN  data  file  is  more  streamlined  than 
the  IDEAS  ASCII  output.  This  both  reduces  the  size  of  the  files  that  have  to  be 
manipulated  and  simplifies  the  post  processing.  A  sample  NASTRAN  output  for  the  six 
tetrahedron  cube  problem  is  shown  in  Section  E.2,  and  the  TetSN  input  file  generated 
from  the  NASTRAN  file  is  shown  in  Section  E.3. 
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E.2  NASTRAN  Output  File 

BEGIN  BULK 


CORD2S* 

2 

0 

0.000000000 

0.000000000+1A 

2 

*1A 

2 

0.000000000 

0.000000000 

0.000000000 

100. 0000000+ IB 

2 

*1B 

2 

100.0000000 

0.000000000 

0.000000000 

+  1C 

? 

*1C 

2 

CCRD2C* 

1 

0 

0.000000000 

0.000000000+1A 

1 

*1A 

1 

0.000000000 

0.000000000 

0.000000000 

100. 0000000+1B 

1 

*IB 

1 

100.0000000 

0.000000000 

0.000000000 

+  1C 

1 

*1C 

1 

GRID* 

1 

0 

-0.500000000 

-0.50000000  0+GA 

1 

*GA 

1 

0.500000000 

0 

GRID* 

2 

0 

-0.500000000 

0 . 500000000+GA 

2 

*GA 

2 

0.500000000 

0 

GRID* 

3 

0 

0.500000000 

-0.50000000  0+GA 

3 

*GA 

3 

0.500000000 

0 

GRID* 

4 

0 

0.500000000 

0. 500000000+GA 

4 

*GA 

4 

0.500000000 

0 

GRID* 

5 

0 

-0.500000000 

-0. 500000000+GA 

5 

*GA 

5 

-0 . 500000000 

0 

GRID* 

6 

0 

-0.500000000 

0. 500000000+GA 

6 

*GA 

6 

-0.500000000 

0 

GRID* 

7 

0 

0.500000000 

-0. 500000000+GA 

7 

*GA 

7 

-0.500000000 

0 

GRID* 

8 

0 

0.500000000 

0. 500000000+GA 

8 

*GA 

8 

-0.500000000 

0 

CTETRA 

1  2 

4  7 

3  5 

CTETRA 

2  2 

4  3 

1  5 

CTETRA 

3  2 

4  1 

2  5 

CTETRA 

4  2 

4  2 

6  5 

CTETRA 

5  2 

4  6 

8  5 

CTETRA 

6  2 

4  8 

7  5 

MAT1* 

1 

2.068000E+9 

0 . 289999992+MA 

1 

*MA 

1 

0.007820000 

0.000011700 

+MB 

1 

*MB 

1 

1.500000E+7 

1 . 500000E+7 

680000.0000 

+MC 

1 

*MC 

1 

PSOLID* 

2 

1 

0 

+  PA 

2 

*  PA 

2 

PARAM 

AUTOS PC  YES 

PARAM 

POST 

-2 

ENDDATA 
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E.3  Sample  TetSN  Input  File 

6  tetrahedron  cube 
# NODES 
8 

X,Y,Z  COORDINATES  OF  EACH  NODE 

-0.500000000000000000  -0.500000000000000000  0.500000000000000000 
-0.500000000000000000  0.500000000000000000  0.500000000000000000 
0.500000000000000000  -0.500000000000000000  0.500000000000000000 
0.500000000000000000  0.500000000000000000  0.500000000000000000 
-0.500000000000000000  -0.500000000000000000  -0.500000000000000000 
-0. 500000000000000000  0. 500000000000000000  -0. 500000000000000000 
0. 500000000000000000  -0.500000000000000000  -0. 500000000000000000 
0.500000000000000000  0.500000000000000000  -0.500000000000000000 
NUMBER  OF  TETRAHEDRA 
6 

NODES  DEFINING  EACH  TETRAHEDRA 

473  5 

43  15 

4  125 

4265 

468  5 

4875 

MATERIAL  NUMBER  OF  EACH  TETRAHEDRA 
1 
1 
1 
1 
1 
1 

NUMBER  OF  MATERIALS 
1 

NUMBER  OF  ENERGY  GROUPS 
1 

SIGMA-TOTAL,  SOURCE,  SIGMA-SCATTER 
1.0000000000000  1 .00000000000000000 
5.0D-1 

BOUNDARY  FACE  DATA 
0 
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